function R=n2rk4(f,a,b,ya,M)
% to use this n2rk4.m you always will need tofind yexact and substitute 

%in this program. ALso you need a file f88.m        which contains three lines
%    function  f=f88(t,y)
%                f=(t-y)/2;
%                end
%
%Input    - f is the function entered as a string 'f'
%         - a and b are the left and right endpoints
%         - ya is the initial condition y(a)
%         - M is the number of steps
%Output  - R = [T' Y' EXA' ERRoR'] where T is the vectof abscissas
%          and Y is the vector of ordinates% modifications to  the code of fink and mathews
% NUMERICAL METHODS: MATLAB Programs
%(c) 1999 by John H. Mathews and Kurtis D. Fink
%To accompany the textbook:
%NUMERICAL METHODS Using MATLAB,
%by John H. Mathews and Kurtis D. Fink
%ISBN 0-13-270042-5, (c) 1999
%PRENTICE HALL, INC.
%Upper Saddle River, NJ 07458
disp('  T    RK4Y   Exact EXact-Rk4       EULERY')
h=(b-a)/M;
T=zeros(1,M+1);
Y=zeros(1,M+1);
EXA=zeros(1,M+1);
ERRO=zeros(1,M+1);
EULR=zeros(1,M+1);
T=a:h:b;
disp('h='); disp(h)
Y(1)=ya;
EULR(1)=ya;
EXA(1)=ya;
ERRO(1)=EXA(1)-Y(1);
for j=1:M
   k1=h*feval(f,T(j),Y(j));
   k2=h*feval(f,T(j)+h/2,Y(j)+k1/2);
   k3=h*feval(f,T(j)+h/2,Y(j)+k2/2);
   k4=h*feval(f,T(j)+h,Y(j)+k3);
   Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)/6;
   EXA(j+1)=T(j+1)-2+3*exp(-1/2*T(j+1));
   ERRO(j+1)=EXA(j+1)-Y(j+1);
   EULR(j+1)=EULR(j)+h*feval(f,T(j),EULR(j));
end

R=[T' Y' EXA' ERRO' EULR'];
T=T';
Y=Y';
clc
xlabel(' T ');
ylabel(' Y ');
title(' RK4-EUlER-EXACT  ');
T=T';
Y=R'
hold
plot(T,Y,'-');
hold
plot(T,EXA,'+')
hold
plot(T,EULR,'.')
return