resol3.m

%rsolution de y'+xy = x

f=inline('[-y(1)./(x+eps)-y(2); y(1)]','x','y');
h = 0.3;
x = 0:h:10;
ye = [0;1];
for i=x(1:end-1)
   ye = [ye euler(f,i,ye(:,end),h)];
end;
y2 = [0;1];
for i=x(1:end-1)
   y2 = [y2 rk2(f,i,y2(:,end),h)];
end;
y4 = [0;1];
for i=x(1:end-1)
   y4 = [y4 rk4(f,i,y4(:,end),h)];
end;
delete(gcf)
plot(x,ye(2,:),'r+')   
hold on
plot(x,y2(2,:),'b*')   
plot(x,y4(2,:),'mo')   
plot(x,besselj(0,x),'y-')
y45 = [0;1];
[x y45]=ode45(f,x,y45);
plot(x,y45(:,2),'g^')
legend('euler','rk2','rk4',...
       'solution analytique','ode45');