[R]=implicit_solver(f,a,ya,b,h,RKa,RKb,RKc,RKns)

function [R]=implicit_solver(f,a,ya,b,h,RKa,RKb,RKc,RKns)
% SOLVER - Returns to the solution of a one-dimensional IVP 
%		   by means of an implicit method choosen by implicit_setmethod. 
%
% CALL [RKa,RKb,RKc,RKns,RKord,RKtype,RKname]=implicit_solver(f,a,ya,b,h,RKa,RKb,RKc,RKns)
%
% INPUT:
%		- f is the function entered as a string 'f'
%       - a and b are the endpoints of the interval
%       - ya is the initial condition y(a)
%       - h is stepsize					
%		- RKa is matrix A from the method's Butcher tableau
%       - RKb is vector b from the method's Butcher tableau
%       - RKc is vector c from the method's Butcher tableau or the
%                 abscissae values in case of Magnus methods
%       - RKns is number of stages in the method or the
%				  number of points method uses
% OUTPUT:
%		- yexact is exact solution vector of ordinates
%		- R is the numerical solution vector of ordinates
%
% Written by Ali Dinler 2006


M=(b-a)/h;
t=zeros(1,M+1);
t=a:h:b;
y(1)=ya;
yexa(1)=ya;
yerr(1)=0;

ini=ones(1,RKns); % initial guess for k1, k2 ,..., ks
jac=zeros(RKns,1);
F=zeros(RKns,1);
max1=100;
epsilon=1e-4;
dh=1e-4;
jacobian=zeros(RKns);
k0=ini';
for i=1:M
   %solve by means of multidimensional newton iteration
   for k=1:max1
    %construct the system    
    for j=1:RKns
       sum1=0;
       for l=1:RKns
           sum1=sum1+RKa(j,l)*k0(l); 
       end
       ff=feval(f,t(i)+RKc(j)*h,y(i)+h*sum1);
       F(j)=k0(j)-ff;
       %obtain jacobian matrix numerically 
       ffh=feval(f,t(i)+RKc(j)*h,y(i)+h*sum1+dh);
       jac(j)=(ffh-ff)/dh;      
       for l=1:RKns
           if(l==j)
            jacobian(j,l)=1-h*RKa(j,l)*jac(j);
           else
            jacobian(j,l)=-h*RKa(j,l)*jac(j);
           end
       end
    end
    temp=jacobian*k0;   
	RHS=-F+temp;
    k1=jacobian\RHS;
    e=k1-k0;
	err=norm(e);
	k0=k1;
    if(err<epsilon),break,end 
   end
   
   sum2=0;
   for l=1:RKns
       sum2=sum2+RKb(l)*k0(l); 
   end 
    
   y(i+1)=y(i)+h*sum2;
   
end

R=y;