f=inline('-2*x*y','x','y');
a=[-1,3];
y0=exp(-1);
n=40;
[y,x]=euler(f,a,y0,n);
figure; plot(x,y);
des=fopen('euler1.txt','w');
fprintf(des,'{');fprintf(des, '{%0.5f, %0.5f},\n', [x; y]); fprintf(des,'}\n');
fclose(des);
[y,x]=modeuler(f,a,y0,n);
figure; plot(x,y);
des=fopen('modeuler1.txt','w');
fprintf(des,'{');fprintf(des, '{%0.5f, %0.5f},\n', [x; y]); fprintf(des,'}\n');
fclose(des);
[y,x]=cranknicholson(f,a,y0,n);
figure; plot(x,y)
des=fopen('cn1.txt','w');
fprintf(des,'{');fprintf(des, '{%0.5f, %0.5f},\n', [x; y]); fprintf(des,'}\n');
fclose(des);
[y,x]=leapfrog(f,a,y0,n);
figure; plot(x,y)
des=fopen('leapfrog1.txt','w');
fprintf(des,'{');fprintf(des, '{%0.5f, %0.5f},\n', [x; y]); fprintf(des,'}\n');
fclose(des);

