
clear all;
a=0.; b=pi;
t0=a; y0=1.;
f=@(t,y) -100*y+100*cos(t)-sin(t);
yd=@(t) cos(t);


method='explicit';
#possbile values for method: 'explicit', 'implicit', 'midpoint', 'trapezoidal'
er=[];
for n=[20,40,80,160,320]
    t=linspace(a,b,n+1);
    h=(b-a)/n;
    
    ydd=[];
    for k =1:n+1
        ydd(k)=yd(t(k));
    end
switch method
  case 'explicit'    
        ya=[];
        ya(1)=y0;
        for k =1:n
            yar=ya(k)+h*f(t(k),ya(k));
            ya=[ya,yar];
        end
    
   case 'implicit'
         ya=[];
        ya(1)=y0;
        for k =1:n
            yar=(ya(k)+h*(100*cos(t(k+1))-sin(t(k+1))))/(1+100*h);%1/(1-h)*ya(k);
            ya=[ya,yar];
        end
            
    case 'trapez'
        ya=[];
        ya(1)=y0;
        for k =1:n
            yar=(1+0.5*h)/(1-0.5*h)*ya(k);
            ya=[ya,yar];
        end
        
    case 'midpoint'
        ya=[];
        ya=[ya,y0];
        ya=[ya,y0+h*f(t0,y0)];
        for k=2:n
            yar=ya(k-1)+2*h*ya(k);
            ya=[ya,yar];
        end
endswitch
    
    figure();
    plot(t,ya,'.',t,ydd,'.');
    
    e=max(abs(ya-ydd));
    er=[er,e];
end
wykl=er(1:length(er)-1)./er(2:length(er));

s="2^a =";
for i=1:length(er)-1
  s=[s," %f,"];
end
s=[s,"\n"];

printf(s, wykl);
