deff('ydot=f(t,y)','ydot=1-2*cos(3*t)*y')
//
function ysol=solalgo(y0,f,T,N)
ysol=[y0];
yn=y0;
for n=0:(N-1)   
deff('yval=eulerimp(y)','yval=y-yn-T/N*f((n+1)*T/N,y)')
yn=fsolve(yn,eulerimp);
ysol=[ysol,yn];
end
endfunction
//
y0=1;T=1;
ysol1=solalgo(y0,f,T,10)
ysol2=solalgo(y0,f,T,100)
ysol3=solalgo(y0,f,T,1000)
yode=ode(y0,0,0:(1/1000):1,f);
figure(1)
clf()
plot2d(0:(1/10):1,ysol1,1);
plot2d(0:(1/100):1,ysol2,2);
plot2d(0:(1/1000):1,ysol3,3);
plot2d(0:(1/1000):1, yode,4);
legends(['N=10','N=100','N=1000','ode'],[1,2,3,4])
