function R=RK42IVPA(g,f,a,b,ya,xa,M);
%- The Initial value problem  y''=f(t,y,y') y(a)=ya, and y'(a)=xa is solvedy Runge Kutta 4 method.
%- Let x' =y = g(t,x,y) then
%-     y'= f(t,x,y) with x(a)=xa, y(a)=ya
%- that means we use RK2
%
% ALso you need two files  f81.m   and f82.m     Each contains three lines
% function       g=f81(t,y,x)
%                g=y
%    function  f=f82(t,y,x)
%                f=y-x^2y-x+0.01*sin(t)      
%                end
% - The exact solution for y''+y=t  with say y(0)=1, and y'(0)=2 is
% -                   yunique=t+cos(t)+sin(t)
%                     xunique=1-sin(t)+cos(t)
%Input    - f is the function entered as a string 'f'
%         -  g   is a a string   'g'
%         - a and b are the left and right endpoints
%         - ya is the initial condition y(a)
%         -yprimea =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

disp('  T    RK4Y RK4X   Exact EXact-Rk4')
h=(b-a)/M;
T=zeros(1,M+1);
Y=zeros(1,M+1);
X=zeros(1,M+1);
EXA=zeros(1,M+1);
EXAP=zeros(1,M+1);
T=a:h:b;
disp('h='); disp(h)
Y(1)=ya;
X(1)=xa
for j=1:M
   k1y=h*feval(f,T(j),Y(j),X(j));
   k1x=h*feval(g,T(j),Y(j),X(j));
   k2y=h*feval(f,T(j)+h/2,Y(j)+k1y/2,X(j)+k1x/2);
   k2x=h*feval(g,T(j)+h/2,Y(j)+k1y/2,X(j)+k1x/2);
   k3y=h*feval(f,T(j)+h/2,Y(j)+k2y/2,X(j)+k2x/2);
   k3x=h*feval(g,T(j)+h/2,Y(j)+k2y/2,X(j)+k2x/2);
   k4y=h*feval(f,T(j)+h,Y(j)+k3y,X(j)+k3x);
   k4x=h*feval(g,T(j)+h,Y(j)+k3y,X(j)+k3x);
   X(j+1)=X(j)+(k1y+2*k2y+2*k3y+k4y)/6;
   Y(j+1)=Y(j)+(k1x+2*k2x+2*k3x+k4x)/6;
   EXA(j+1)=T(j+1)+cos(T(j+1))+sin(T(j+1));
   EXAP(j+1)=1-sin(T(j+1))+cos(T(j+1));
end
R=[T' Y' X' EXA' EXAP'];
T=T';
X=X';
Y=Y';
EXA=T'+cos(T')+sin(T');
EXAP=1+sin(T')+cos(T');
clc
xlabel(' T ');
ylabel(' Y ');
title(' RK4Y,X -IVP-EXACT  ');
hold on
plot(T,Y,'>');
plot(T,EXA,'+');
plot(T,EXAP,'-');
plot(T,X,'<');
return