function [y_approx,Err,w,t,m]=RK4(f,a,b,y_in,y_ex,N,plott,error)
%------------------------------------------------
%function RK4:
%
%   INPUT:  a     the end point
%           b     the other end point
%           N     max number of steps
%           TOL   tolerance
%           y_in  initial value at y(a)
%           y_ex  exact value at y(b)
%           plott if 1 is entered, the program will plot
%           error if 1 is entered, the program will calculate error based on y_ex
%           f     the input functions
%
%   OUTPUT: y_ap  approximate answer for y(b)
%           w     computed solution
%           Err   error between the exact value of
%                 our calculated value
%           m     number of equations
%           t     time array
%
%
% example: mid_pt('f',0,2,[1 0 0],[25 25 25],1000,1,1); where f is f.m
%-------------------------------------------------

t(1)=a;
w(1,:)=y_in;                % set w(0) = the initial value
m=length(y_in);             % set the number of equations
h=(b-a)/N;                  % set h

for i=1:N
    k1 = h*feval(f, t(i), w(i,:));
    k2 = h*feval(f, (t(i)+h/2), (w(i,:)+k1/2));
    k3 = h*feval(f, (t(i)+h/2), (w(i,:)+k2/2));
    k4 = h*feval(f, (t(i)+h), (w(i,:)+k3));
    w(i+1,:) = w(i,:)+(k1+2*k2+2*k3+k4)/6;
    t(i+1) = t(i)+h;
end
y_approx=w(N+1,:);
Err=abs(y_ex - y_approx);

if plott==1
    figure
    hold
    for j=1:m
        plot(t,w(:,j))
        
    end;
    
end;

if error==1
    % output stuff...
    fprintf('\n');
    disp('mid - point method');
    disp('------------------------------');
    disp(sprintf('exact Val:\t%f',y_ex));
    disp(sprintf('Cal. Val:\t%f',y_approx));
    disp(sprintf('Error:  \t%d',Err));
    disp(sprintf('Step:   \t%d',N));
    fprintf('\n');
end