function numsolver(action)
% NUMSOLVER(action) is a graphical tool to explore 2nd order Differential 
%           Equations.  In order to use the tool, a 2nd order ODE must be written
%		 	   in terms of a 2nd order system.  Currently, both Euler's method and 
%			   4th order Runga Kutta method's are available.

 global tnum xnum ynum;
 
%   if isempty(action)
%     numsolver('start');
%   end


% These few lines let's the user just type "numsolver" to run the program.

   if nargin ==0 
      numsolver('start');
   else
   switch action	
      
      %  This code reads in all of the data enter through the window and runs the
      %  numerical solver.
      
	case 'runsolver'
      
      % This is the first and second equation, stored as a String
      
    xEqn=get(findobj('Tag','xEqnTextBox'),'String');
    yEqn=get(findobj('Tag','yEqnTextBox'),'String');
    
    % Open the file 'oderhs.m' to write the equation to it.
    
    fid=fopen('oderhs.m','w');
    
    % Create the code for the function that will be called by the solver.
    
    fprintf(fid,'function yprime=oderhs(time,variable)\n');
    fprintf(fid,'\n\n x=variable(1);y=variable(2);t=time;\n');
    fprintf(fid,['\n\n yprime=[' xEqn ';' yEqn '];']);
    fclose(fid);
    
    % t0 is the initial time
    % dt is the timestep
    % x0 is the initial value for x (the first variable)
    % y0 is the initial value for y (the second variable)
    % numsteps is the number of steps of the solver to use.
    
    t0=str2num(get(findobj('Tag','t0TextBox'),'String'));
    dt=str2num(get(findobj('Tag','StepSizeTextBox'),'String'));
    x0=str2num(get(findobj('Tag','x0TextBox'),'String'));
    y0=str2num(get(findobj('Tag','y0TextBox'),'String'));
    numsteps=str2num(get(findobj('Tag','NumStepsTextBox'),'String'));
    
    % The variable SCMenu is where the information about while solver is stored
    %  if get(SCMenu,'Value')==1 then 4th order Runga Kutta
    %  if get(SCMenu,'Value')==2 then Euler's Method.
    
    SCMenu=findobj('tag','SolverChoice');
    
    if get(SCMenu,'Value')==1
       HF_wait=waitbar(0,'Solving ODE using Runga Kutta.  Please Wait');   
    else
       HF_wait=waitbar(0,'Solving ODE using Euler''s.  Please Wait');   
    end
    
    x(1)=x0;y(1)=y0;t(1)=t0;
    
    if get(SCMenu,'Value')==1  % The Runga Kutta Method
        for i=1:numsteps
         
         waitbar(i/numsteps);
         
         t(i+1)=t0+i*dt;
         
         % The following lines of code is the actual RK4 method.
         
         out=dt*feval('oderhs',t(i),[x(i);y(i)]);
         x1=out(1);y1=out(2);
         out=dt*feval('oderhs',t(i)+0.5*dt,[x(i)+0.5*x1;y(i)+0.5*y1]);
         x2=out(1);y2=out(2);	
         out=dt*feval('oderhs',t(i)+0.5*dt,[x(i)+0.5*x2;y(i)+0.5*y2]);
         x3=out(1);y3=out(2);
			out=dt*feval('oderhs',t(i)+dt,[x(i)+x3;y(i)+y3]);
         x4=out(1);y4=out(2);
         
 	 		x(i+1)=x(i)+(x1+2*x2+2*x3+x4)/6;
	 		y(i+1)=y(i)+(y1+2*y2+2*y3+y4)/6;	 
      end
      
      
   else  % use the Euler's method
      	t(1)=t0;
      for i=1:numsteps 
         waitbar(i/numsteps);
         
         t(i+1)=t0+dt*i;
			out=dt*feval('oderhs',t(i),[x(i);y(i)]);
         
         % Here's Euler's Method
         
         x(i+1)=x(i)+out(1);
         y(i+1)=y(i)+out(2);
         
      end
   end
   
   close(HF_wait);
   
   TSFig=findobj('tag','TSFig');  % The Handle for the Time Series Figure
    
    if isempty(TSFig)   % If the window doesn't exist create it.
      TSFig=figure('tag','TSFig','Position',[10 400 400 300]);
      hold on	             % This ensures that subsequent plots are added 
       	 						 % to the axes.  
      plot(t,x,'g-',t,y,'b-.');  %Plot the results
      legend('x(t)','y(t)');     % Add a legend to the plot      
      xlabel('t');ylabel('x(t), y(t)');  % label the axes
    else
       figure(TSFig);
       plot(t,x,'g-',t,y,'b-.');  % Plot the results
       legend('x(t)','y(t)');	    % Add a legend
       xlabel('t'); ylabel('x(t), y(t)');  % label the axes
    end
    
    PPFig=findobj('tag','PPFig');  % The handle for the Phase Portrait
    
    if isempty(PPFig)  % If the window doesn't exist, create it.
       PPFig=figure('tag','PPFig', 'Position',[10 50 400 300]);
       plot(x,y);           % Plot the two variables with respect to each other
       hold on;				 % This ensures that subsequent plots are added 
       	 						 % to the axes.  
       xlabel('x(t)');ylabel('y(t)');
    else
       figure(PPFig);
       plot(x,y);
       xlabel('x(t)');ylabel('y(t)');
    end
    
    clear oderhs;  % This clear the function oderhs from the memory, so if the
      			    % the user changes the ODE, then it must be rebuilt
    
    % end the runsolver part of the function.

    case 'clearplots'   % clears the plots
    
     TSFig=findobj('tag','TSFig'); 
     PPFig=findobj('tag','PPFig'); 
     
     
     if ~isempty(TSFig)
       figure(TSFig)
       cla;
     end
     
     if ~isempty(PPFig)
        figure(PPFig)   
	cla;
     end
      
   
   % end the clearplots part of the function 
    
  case 'closesolver'  % close the numsolver GUI.
    
     SolverFig=findobj('Tag','EulerParamFig');
     TSFig=findobj('tag','TSFig'); 
     PPFig=findobj('tag','PPFig'); 
     close(SolverFig);
     
     
     if ~isempty(TSFig)
        close(TSFig)
     end
     
     if ~isempty(PPFig)
        close(PPFig)   
     end
     
     % end the closesolver part of the function
     
   case 'start'  % The code to create the numsolver GUI.
      
      % caution:  This code is long and boring.  It was created using the "guide" 
      % function in Matlab.  Proceed carefully.
      
      
    h0 = figure('Units','points', ...
	'Color',[0.8 0.8 0.8], ...
	'Name','Numerical Solver for 2D systems', ...
	'NumberTitle','off', ...
	'Position',[435.75 251.25 322.5 269.25], ...
	'Tag','EulerParamFig');
    h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.8 0.8 0.8], ...
	'FontSize',12, ...
	'ListboxTop',0, ...
	'Position',[31.5 188.25 21.75 16.5], ...
	'String','dy', ...
	'Style','text', ...
	'Tag','StaticText1');
    h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.8 0.8 0.8], ...
	'FontSize',12, ...
	'ListboxTop',0, ...
	'Position',[30.75 222 21.75 17.25], ...
	'String','dx', ...
	'Style','text', ...
	'Tag','StaticText1');
    h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.8 0.8 0.8], ...
	'FontSize',12, ...
	'ListboxTop',0, ...
	'Position',[30.9512 206.674 21.9654 17.9717], ...
	'String','dt', ...
	'Style','text', ...
	'Tag','StaticText1');
    h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0 0 0], ...
	'ListboxTop',0, ...
	'Position',[29.25 225.75 22.5 2.25], ...
	'Style','frame', ...
	'Tag','Frame1');
    h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.8 0.8 0.8], ...
	'FontSize',14, ...
	'FontWeight','bold', ...
	'ListboxTop',0, ...
	'Position',[57 216.75 22.5 19.5], ...
	'String','=', ...
	'Style','text', ...
	'Tag','StaticText1');
    h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[1 1 1], ...
	'FontSize',12, ...
	'HorizontalAlignment','left', ...
	'ListboxTop',0, ...
	'Position',[74.25 219 220.5 21], ...
	'String','-y', ...
	'Style','edit', ...
	'Tag','xEqnTextBox');
    h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.8 0.8 0.8], ...
	'FontSize',14, ...
	'FontWeight','bold', ...
	'ListboxTop',0, ...
	'Position',[55.5 177 22.5 18.75], ...
	'String','=', ...
	'Style','text', ...
	'Tag','StaticText1');
    h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0 0 0], ...
	'ListboxTop',0, ...
	'Position',[28.5 186 24.75 2.25], ...
	'Style','frame', ...
	'Tag','Frame1');
    h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.8 0.8 0.8], ...
	'FontSize',12, ...
	'ListboxTop',0, ...
	'Position',[28.5 174.75 24.75 10.5], ...
	'String','dt', ...
	'Style','text', ...
	'Tag','StaticText1');
    h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[1 1 1], ...
	'FontSize',12, ...
	'HorizontalAlignment','left', ...
	'ListboxTop',0, ...
	'Position',[72 178.5 221.25 20.25], ...
	'String','x', ...
	'Style','edit', ...
	'Tag','yEqnTextBox');
    h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.8 0.8 0.8], ...
	'FontSize',12, ...
	'HorizontalAlignment','left', ...
	'ListboxTop',0, ...
	'Position',[19.5 121.5 77.25 19.5], ...
	'String','x(t0)', ...
	'Style','text', ...
	'Tag','StaticText1');
    h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.8 0.8 0.8], ...
	'FontSize',12, ...
	'HorizontalAlignment','left', ...
	'ListboxTop',0, ...
	'Position',[20.9669 90.8567 76.8787 19.9685], ...
	'String','y(t0)', ...
	'Style','text', ...
	'Tag','StaticText1');
    h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.8 0.8 0.8], ...
	'FontSize',12, ...
	'HorizontalAlignment','left', ...
	'ListboxTop',0, ...
	'Position',[93 123 77.25 19.5], ...
	'String','Initial time (t0)', ...
	'Style','text', ...
	'Tag','StaticText1');
    h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.8 0.8 0.8], ...
	'FontSize',12, ...
	'HorizontalAlignment','left', ...
	'ListboxTop',0, ...
	'Position',[93 93.75 61.5 19.5], ...
	'String','Step Size', ...
	'Style','text', ...
	'Tag','StaticText1');
    h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.8 0.8 0.8], ...
	'FontSize',12, ...
	'HorizontalAlignment','left', ...
	'ListboxTop',0, ...
	'Position',[91.5 68.25 60.75 19.5], ...
	'String','Num. Steps', ...
	'Style','text', ...
	'Tag','StaticText1');
    h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[1 1 1], ...
	'FontSize',12, ...
	'HorizontalAlignment','left', ...
	'ListboxTop',0, ...
	'Position',[44.25 122.25 44.25 21], ...
	'String','1', ...
	'Style','edit', ...
	'Tag','x0TextBox');
    h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[1 1 1], ...
	'FontSize',12, ...
	'HorizontalAlignment','left', ...
	'ListboxTop',0, ...
	'Position',[41.25 93.75 47.25 20.25], ...
	'String','0', ...
	'Style','edit', ...
	'Tag','y0TextBox');
    h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[1 1 1], ...
	'FontSize',12, ...
	'HorizontalAlignment','left', ...
	'ListboxTop',0, ...
	'Position',[156 127.5 49.5 15], ...
	'String','0', ...
	'Style','edit', ...
	'Tag','t0TextBox');
    h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[1 1 1], ...
	'FontSize',12, ...
	'HorizontalAlignment','left', ...
	'ListboxTop',0, ...
	'Position',[156 100.5 49.5 15], ...
	'String','0.1', ...
	'Style','edit', ...
	'Tag','StepSizeTextBox');
    h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[1 1 1], ...
	'FontSize',12, ...
	'HorizontalAlignment','left', ...
	'ListboxTop',0, ...
	'Position',[156 75 49.5 15], ...
	'String','100', ...
	'Style','edit', ...
	'Tag','NumStepsTextBox');
    h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'Callback','numsolver(''runsolver'');', ...
	'FontSize',12, ...
	'FontWeight','bold', ...
	'ListboxTop',0, ...
	'Position',[10.75 33.75 83.25 32.25], ...
	'String','Solve', ...
	'Tag','SolveButton');
    h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'Callback','numsolver(''clearplots'');', ...
	'FontSize',12, ...
	'FontWeight','bold', ...
	'ListboxTop',0, ...
	'Position',[110.75 33.75 83.25 32.25], ...
	'String','Clear', ...
	'Tag','CloseButton'); 
    h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'Callback','numsolver(''closesolver'');', ...
	'FontSize',12, ...
	'FontWeight','bold', ...
	'ListboxTop',0, ...
	'Position',[212.75 33.75 83.25 32.25], ...
	'String','Close', ...
	'Tag','CloseButton');
    h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'ListboxTop',0, ...
	'Position',[226.5 99.75 75 17.25], ...
	'String',[' Runga Kutta ';' Euler''s     '], ...
	'Style','popupmenu', ...
	'Tag','SolverChoice', ...
	'Value',1);
    h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.8 0.8 0.8], ...
	'FontSize',12, ...
	'HorizontalAlignment','left', ...
	'ListboxTop',0, ...
	'Position',[222 124.5 77.25 19.5], ...
	'String','Solver Type', ...
	'Style','text', ...
	'Tag','StaticText1');
    if nargout > 0, fig = h0; end

end   % ends the switch-case flow

end	% ends the if nargin==0 flow.