function [soln,error,order]=numsolve(FFFCCNNN,WINvect,derivstr,noSteps,Tname,fieldLine)
% NUMSOLVE  solving equations numerically and plots the solution.
%
% usage:
%   [soln,error,order] = NUMSOLVE(FFFCCNNN, WINvect, derivstr, noSteps, Tname, fieldLine)
%   FFFCCNNN = name of the function file describing the problem
%   WINvect = 1x5 vector containing following informations -
%             WINvect(1) = initial value of independent variable
%             WINvect(2) = final value of independent variable
%             WINvect(3) = value of dependent variable at initial independent variable
%             WINvect(4) = (if known) value of dependent variable at final independent variable,
%                          if the value is unidentified, Inf can be used instead
%             WINvect(5) = method to be used for the calculation
%                        = 1   Euler
%                        = 2   Improved Euler
%                        = 3   Runge-Kutta
%   derivstr = string of differential equation described by FFFCCNNN function without
%              'dx/dt' bit. eg. if dx/dt = 3*t is to be calculated, derivstr = '3*t'
%   noSteps = number of steps to be calculated to. (in powers of 2)
%   Tname = string specifies the name of the independent variable (eg. 't')
%   fieldLine = 1x2 vector contains information about the field lines on the plot
%               fieldLine(1) = whether to draw field lines, 0 if no, 1 if yes
%               fieldLine(2) = number of line segments to be drawn
%
% Modified and updated version of written program.
% Modified and updated by Il-Joong Kil, Febuary 2002
%
% Bug fixed by Mike Meylan March 2005.  Please report any bugs to 
% meylan@math.auckland.ac.nz.  


steps = noSteps;
m     = WINvect(5);
order = zeros(1,10);
error = zeros(1,10);
tdraw = 0; ydraw = 0;
if m==1
    meth='Euler';
    for i=1:10,
        t=[];y=[];
        [t y]=eul(FFFCCNNN,[WINvect(1) WINvect(2)],WINvect(3),(-WINvect(1) + WINvect(2))/(2^(i-1)));
        if(i == steps + 1)
            tdraw = t; ydraw = y;
        end
        soln(i)=y(length(t));
    end
elseif m==2
    for i=1:10,
        meth='Improved Euler';
        t=[];y=[];
        [t y]=rk2(FFFCCNNN,[WINvect(1) WINvect(2)],WINvect(3),(WINvect(1) - WINvect(2))/(2^(i-1)));
        if(i == steps + 1)
            tdraw = t; ydraw = y;
        end
        soln(i)=y(length(t));
    end  
else
    for i=1:10
        meth='4th order Runge-Kutta';
        t=[];y=[];
        [t y]=rk4(FFFCCNNN,[WINvect(1) WINvect(2)],WINvect(3),(WINvect(1) - WINvect(2))/(2^(i-1)));
        if(i == steps + 1)
            tdraw = t; ydraw = y;
        end
        soln(i)=y(length(t));
    end
end;
t = tdraw; y = ydraw;
new_fig = figure;
set(new_fig,'name','NUMERICAL PLOT','numb','off','menubar','none')
set(new_fig,'CloseRequestFcn',['warndlg(cellstr({''You can only close this figure by clicking ''''Close'''' button'',',...
        '''from the corresponding NUMERICAL Display window or'',',...
        '''from main interface by clicking ''''Close Plots'''' button''}))']);
f = uimenu('Label','&Menu');
uimenu(f,'Label','&Copy to clipboard','Callback','print -dbitmap',...
    'Accelerator','C');
uimenu(f,'Label','&Print...','Callback','printdlg','Accelerator','P');
uimenu('Label','&Grid on/off','Callback','grid');
g = uimenu('Label','&Zoom');
zoomCmd1 = ['h = findobj(gcf,''Label'',''Zoom &in'');',...
        'if strcmp(get(h,''checked''),''off'') ','zoom on;',...
        'set(h,''Checked'',''on'');','else zoom off;',...
        'set(h,''Checked'',''off''); end;'];
zoomCmd2 = ['zoom off;' 'zoom out;' 'h = findobj(gcf,''Label'',''Zoom &in'');',...
        'set(h,''Checked'',''off'');'];
uimenu(g,'Label','Zoom &in','checked','off','callback',zoomCmd1);
uimenu(g,'Label','Zoom to &original size','callback',zoomCmd2);
plot(t,y); hold on; ax = axis;
if (fieldLine(1) == 1) plotdefield(FFFCCNNN,ax,fieldLine(2)); end
axis(ax);
xlabel(Tname); ylabel('x');
h=abs(WINvect(2)-WINvect(1))/(2^(steps));stepsize=['h = ',num2str(h)];
if length(Tname)>1 hstr = ['dx/d(' Tname ')'];
else hstr = ['dx/d' Tname]; end;
temp1=['Numerical solution to \bf',hstr,' = ',derivstr,'\rm drawn after performing ' num2str(2^steps) ' steps'];
temp2=['\bf' meth, '\rm method is used, with stepsize \bf' stepsize];
t1 = length(temp1); t2 = length(temp2); td = abs(t1-t2);
if t1 < t2
    if rem(td,2) == 0
        for i = 1:td/2 temp1 = [' ' temp1 ' ']; end;
    else temp1 = [temp1 ' '];
        for i = 1:(td-1)/2 temp1 = [' ' temp1 ' ']; end;
    end
elseif t1 > t2
    if rem(td,2) == 0
        for i = 1:td/2 temp2 = [' ' temp2 ' ']; end;
    else temp2 = [temp2 ' '];
        for i = 1:(td-1)/2 temp2 = [' ' temp2 ' ']; end;
    end
end
header = [temp1;temp2]; ax = axis;
title(header,'Position',[(ax(2)+ax(1))/2 ax(4) 0])

if WINvect(4)==inf,
    error = inf*ones(1,10);
    for k=1:8,
        order(k)=log((soln(k+1)-soln(k))/(soln(k+2)-soln(k+1)))/log(2);
    end;
else
    for i=1:10,
        error(i)=abs(soln(i)-WINvect(4));
    end;
    for i=1:9,
        order(i)=(log(error(i))-log(error(i+1)))/log(2);
    end; 
    
end;

for k = 1:length(order)
    % 4 different cases for the effective order(EO) to be undefined
    if isreal(order(k)) & real(order(k))~=order(k) order(k) = 1*j;  % EO is 0/0
    elseif order(k)==Inf order(k) = 2*j;                            % EO is a/0, a>0
    elseif order(k)==-Inf order(k) = 3*j;                           % EO is a/0, a<0
    elseif real(order(k))~=order(k) order(k) = 4*j;                 % EO is in complex form*
    end
end
% * due to the nature of log function in MatLab which gives
% complex number as the answer to log of negative number.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [tout, yout] = eul(FunFcn, tspan, y0, ssize)

% EUL	Integrates a system of ordinary differential equations using
%	Euler's method.  See also ODE45 and ODEDEMO.
%	[t,y] = eul('yprime', tspan, y0) integrates the system
%	of ordinary differential equations described by the M-file
%	yprime.m over the interval tspan = [t0,tfinal] and using initial
%	conditions y0.
%	[t, y] = eul(F, tspan, y0, ssize) uses step size ssize
%
% INPUT:
% F     - String containing name of user-supplied problem description.
%         Call: yprime = fun(t,y) where F = 'fun'.
%         t      - Time (scalar).
%         y      - Solution vector.
%         yprime - Returned derivative vector; yprime(i) = dy(i)/dt.
% tspan = [t0, tfinal], where t0 is the initial value of t, and tfinal is
%         the final value of t.
% y0    - Initial value vector.
% ssize - The step size to be used. (Default: ssize = (tfinal - t0)/100).
%
% OUTPUT:
% t  - Returned integration time points (column-vector).
% y  - Returned solution, one solution row-vector per tout-value.
%
% The result can be displayed by: plot(t,y).

% Initialization

t0=tspan(1);
tfinal=tspan(2);
pm = sign(tfinal - t0);  % Which way are we computing?
if (nargin < 4), ssize = abs(tfinal - t0)/100; end
if ssize < 0, ssize = -ssize; end
h = pm*ssize;  
t = t0;
y = y0(:);
tout = t;
yout = y.';

% We need to compute the number of steps.

dt = abs(tfinal - t0);
N = floor(dt/ssize) + 1;
if (N-1)*ssize < dt
  N = N + 1;
end

% Initialize the output.

tout = zeros(N,1);
tout(1) = t;
yout = zeros(N,size(y,1));
yout(1,:) = y.';
k = 1;

% The main loop

while k < N
  if pm*(t + h - tfinal) > 0
    h = tfinal - t;
    tout(k+1) = tfinal;
  else
    tout(k+1) = t0 +k*h;
  end
  k = k+1;
  % Compute the slope
  s1 = feval(FunFcn, t, y); s1 = s1(:); % s1=f(t(k),y(k))  
  y = y + h*s1;		% y(k+1) = y(k) + h*f(t(k),y(k))
  t = tout(k);
  yout(k,:) = y.';
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [tout, yout] = rk2(FunFcn, tspan, y0, ssize)

% RK2	Integrates a system of ordinary differential equations using
%	the second order Runge-Kutta  method.  See also ODE45 and
%	ODEDEMO.M.
%	[t,y] = rk2('yprime', tspan, y0) integrates the system
%	of ordinary differential equations described by the M-file
%	yprime.m over the interval tspan=[t0,tfinal] and using initial
%	conditions Y0.
%	[t, y] = rk2(F, tspan, y0, ssize) uses step size ssize
%
% INPUT:
% F     - String containing name of user-supplied problem description.
%         Call: yprime = fun(t,y) where F = 'fun'.
%         t      - Time (scalar).
%         y      - Solution column-vector.
%         yprime - Returned derivative column-vector; yprime(i) = dy(i)/dt.
% tspan = [t0, tfinal], where t0 is the initial value of t, and tfinal is
%         the final value of t.
% y0    - Initial value column-vector.
% ssize - The step size to be used. (Default: ssize = (tfinal - t0)/100).
%
% OUTPUT:
% t  - Returned integration time points (column-vector).
% y  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(t,y).


% Initialization

t0=tspan(1);
tfinal=tspan(2);
pm = sign(tfinal - t0);  % Which way are we computing?
if nargin < 4, ssize = (tfinal - t0)/100; end
if ssize < 0, ssize = -ssize; end
h = pm*ssize;
t = t0;
y = y0(:);

% We need to compute the number of steps.

dt = abs(tfinal - t0);
N = floor(dt/ssize) + 1;
if (N-1)*ssize < dt
  N = N + 1;
end

% Initialize the output.

tout = zeros(N,1);
tout(1) = t;
yout = zeros(N,size(y,1));
yout(1,:) = y.';
k = 1;

% The main loop

while k < N
  if pm*(t + h - tfinal) > 0 
    h = tfinal - t; 
    tout(k+1) = tfinal;
  else
    tout(k+1) = t0 +k*h;
  end
  k = k + 1;
  % Compute the slopes
  s1 = feval(FunFcn, t, y); s1 = s1(:);
  s2 = feval(FunFcn, t + h, y + h*s1); s2=s2(:);
  y = y + h*(s1 + s2)/2;
  t = tout(k);
  yout(k,:) = y.';
end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [tout, yout] = rk4(FunFcn, tspan, y0, ssize)

% RK4	Integrates a system of ordinary differential equations using
%	the fourth order Runge-Kutta  method.  See also ODE45 and
%	ODEDEMO.M.
%	[t,y] = rk4('yprime', tspan, y0) integrates the system
%	of ordinary differential equations described by the M-file
%	yprime.m over the interval tspan=[t0,tfinal] and using initial
%	conditions y0.
%	[t, y] = rk4(F, tspan, y0, ssize) uses step size ssize
%
% INPUT:
% F     - String containing name of user-supplied problem description.
%         Call: yprime = fun(t,y) where F = 'fun'.
%         t      - Time (scalar).
%         y      - Solution column-vector.
%         yprime - Returned derivative column-vector; yprime(i) = dy(i)/dt.
% tspan = [t0, tfinal], where t0 is the initial value of t, and tfinal is
%         the final value of t.
% y0    - Initial value column-vector.
% ssize - The step size to be used. (Default: ssize = (tfinal - t0)/100).
%
% OUTPUT:
% t  - Returned integration time points (column-vector).
% y  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(t,y).


% Initialization

t0=tspan(1);
tfinal=tspan(2);
pm = sign(tfinal - t0);  % Which way are we computing?
if nargin < 4, ssize = (tfinal - t0)/100; end
if ssize < 0, ssize = -ssize; end
h = pm*ssize;
t = t0;
y = y0(:);

% We need to compute the number of steps.

dt = abs(tfinal - t0);
N = floor(dt/ssize) + 1;
if (N-1)*ssize < dt
  N = N + 1;
end

% Initialize the output.

tout = zeros(N,1);
tout(1) = t;
yout = zeros(N,size(y,1));
yout(1,:) = y.';
k = 1;

% The main loop
while (k < N)
  if pm*(t + h - tfinal) > 0 
    h = tfinal - t; 
    tout(k+1) = tfinal;
  else
    tout(k+1) = t0 +k*h;
  end
  k = k + 1;

  % Compute the slopes
  s1 = feval(FunFcn, t, y); s1 = s1(:);
  s2 = feval(FunFcn, t + h/2, y + h*s1/2); s2=s2(:);
  s3 = feval(FunFcn, t + h/2, y + h*s2/2); s3=s3(:);
  s4 = feval(FunFcn, t + h, y + h*s3); s4=s4(:);
  y = y + h*(s1 + 2*s2 + 2*s3 +s4)/6;
  t = tout(k);
  yout(k,:) = y.';
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%