function output = numerical(workingDir,action,input1)
% Numerical Approximation.
%   Performs the numerical approximation by Euler, Improved Euler
%   and 4th order Runge-Kutta methods.
%
%   Note that the print button on the solution output only works
%   for the system with less than 32 bit color depth.
%   It will not print properly if the color depth is 32 bit
%
%   Original programmatic layout is adapted from dfield written by John C. Polking, Rice University
%
%   Modified and updated version of written program.
%   Modified and updated by Il-Joong Kil, November 2001
%   Finalised to version 1.7 on 13 March 2002
%
%   Bug fixed in numsolve by Mike Meylan, Marc 2005. Please report any
%   bugs to meylan@math.auckland.ac.nz

format long;

% Routine to decide the working directory for the program
defdir = cd;
if nargin == 0 | nargin == 1
    action ='initialize';
    if nargin == 0 workingDir = defdir; end;
end
try
    cd(workingDir);
catch
    try 
        workingDir = defdir;
        cd(workingDir);
    catch
        workingDir = cd;
        cd(workingDir);
    end
end

% Initialisation sequences
if strcmpi(action,'initialize')
    
    disp(['Directory in use : ', workingDir]);
    figs = get(0,'children');
    for (figno = 1:length(figs))
        nn = get(figs(figno),'name');
        TTT = findstr(nn,'NUMERICAL');
        SSS = findstr(nn,'WARNING!!');
        if (~isempty(TTT)|~isempty(SSS)) delete(figs(figno)); end
    end
    
    % Initiate several parameters.
    ysoln    = zeros(1,10);
    error    = zeros(1,10);
    order    = zeros(1,10);
    
    % The default differential equation, and Display window.
    Tname    = 't'; derivstr = 'x^2 - t'; WINvect  = [0 1 0 inf 1];
    
    dfset = figure('pos',[30 70 470 414],'resize','off','name','NUMERICAL approximation',...
        'numb','off','visible','off','color',[0.95 0.85 0.75],'menu','none');
    figure(dfset)
    
    frame(1) = uicontrol('style','frame','backgroundcolor',[0.8 0.8 0.9],...
        'pos',[25 305 420 90.5],'visible','off');
    
    tname = [
        'sh = get(gcf,''user'');','sm = get(sh(1),''user'');'...
            'eq = sh(4:9);','wind =sh(10:18);','method = sh(22:24);'...
            'Tname=deblank2(get(eq(6),''string''));'...
            'set(eq(6),''string'',Tname);'...
            'if(length(Tname) > 1) set(eq(2),''string'',[''dx/d('', Tname,'')'']);'...
            'else set(eq(2),''string'',[''dx/d'', Tname]); end;'...
            'set(wind(2),''string'',[''initial value of '',Tname,'' = '']);'...
            'set(wind(6),''string'',[''end value of '',Tname,'' = '']);'...
            'set(wind(8),''string'',[''solution at final '',Tname,'' = '']);'...
            'sm = str2mat(Tname,sm(2,:));'...
            'set(sh(1),''user'',sm);'];
    
    xder =[
        'sh = get(gcf,''user'');','sm = get(sh(1),''user'');','eq=sh(4:9);'...
            'derivstr = deblank2(get(eq(4),''string''));'...
            'sm = str2mat(sm(1,:),derivstr);','set(eq(4),''string'',derivstr);'...
            'set(sh(1),''user'',sm);'];
    
    eq(1)=uicontrol('style','text','pos',[30 368 410 20],'fontweight','bold','horizon','center',...
        'string','Enter the differential equation:','visible','off','backgroundcolor',[0.8 0.8 0.9]);
    
    eq(2) = uicontrol('style','text','pos',[30 342 80 20],'horizon','right',...
        'string',['dx/d' Tname],'visible','off','backgroundcolor',[0.8 0.8 0.9]);
    
    eq(3) = uicontrol('style','text','pos',[110 342 30 20],'horizon','center',...
        'string','= ','visible','off','backgroundcolor',[0.8 0.8 0.9]);
    
    eq(4)=uicontrol('pos',[140 345 270 20],'backgroundcolor','w','string',derivstr,...
        'horizon','left','style','edit','call',xder,'visible','off');
    
    eq(5)=uicontrol('pos',[30 307.5 190 20],'style','text','horizon','right',...
        'string','The independent variable is ','visible','off','backgroundcolor',[0.8 0.8 0.9]);
    
    eq(6)=uicontrol('pos',[230 310.5 80 20],'string',Tname,'backgroundcolor','w',...
        'horizon','center','style','edit','call',tname,'visible','off');
    
    for i = 1:3
        w(i,:) = [
            'sh = get(gcf,''user'');','wind = sh(10:18);'...
                'text = deblank2(get(wind(' num2str(i*2+1) '),''string''));'...
                'sh(' num2str(20+i) ')=str2num(text);','set(wind(' num2str(i*2+1) '),''string'',text);'...
                'set(gcf,''user'',sh);'];
    end
    w4 = [
        'sh = get(gcf,''user'');','wind = sh(10:18);'...
            'text = get(wind(' num2str(9) '),''string'');'...
            'if (isempty(deblank(text))) text = ''Inf''; end;'...
            'text = deblank2(text);'...
            'sh(' num2str(24) ')=str2num(text);','set(wind(' num2str(9) '),''string'',text);'...
            'set(gcf,''user'',sh);'];
    
    
    frame(2) = uicontrol('style','frame','pos',[25 195 420 97],'visible','off','backgroundcolor',[0.8 0.8 0.9]);
    
    wind(1)=uicontrol('style','text','fontweight','bold','backgroundcolor',[0.8 0.8 0.9],'horizon','center',...
        'pos',[30 263 410 20],'string','Enter the initial and final conditions','visible','off');
    
    wind(2) = uicontrol('style','text','backgroundcolor',[0.8 0.8 0.9],'pos',[30 232 125 20],...
        'horizon','right','string',['initial value of ',Tname,' = '],'visible','off');
    
    wind(3) = uicontrol('style','edit','backgroundcolor','w','pos',[155 234 50 20],...
        'string',num2str(WINvect(1)),'call',w(1,:),'visible','off');
    
    wind(4) = uicontrol('style','text','backgroundcolor',[0.8 0.8 0.9],'pos',[245 232 115 20],...
        'horizon','right','string',['initial value of x = '],'visible','off');
    
    wind(5)=uicontrol('style','edit','backgroundcolor','w','pos',[155 211 50 20],...
        'string',num2str(WINvect(2)),'call',w(2,:),'visible','off');
    
    wind(6) = uicontrol('style','text','backgroundcolor',[0.8 0.8 0.9],'pos',[30 207 125 20],...
        'horizon','right','string',['end value of ',Tname,' = '],'visible','off');
    
    wind(7) = uicontrol('style','edit','backgroundcolor','w','pos',[360 234 50 20],...
        'string',num2str(WINvect(3)),'call',w(3,:),'visible','off');
    
    wind(8) = uicontrol('style','text','backgroundcolor',[0.8 0.8 0.9],'pos',[225 207 135 20],...
        'horizon','right','string',['solution at final ',Tname,' = '],'visible','off');
    
    wind(9)=uicontrol('style','edit','backgroundcolor','w','pos',[360 211 70 20],...
        'string',num2str(WINvect(4)),'call',w4,'visible','off');
    
    frame(3) = uicontrol('style','frame','pos',[25 63 420 120],'visible','off','backgroundcolor',[0.8 0.8 0.9]);
    
    pop(1) = uicontrol('Style','Popup','pos',[235 150 140 20],'String','Euler|Improved Euler|4th order Runge-Kutta',...
        'Callback',['sh = get(gcf,''user'');' 'sh(25) = get(sh(58),''value'');'],'backgroundcolor','w','value',WINvect(5));
    
    pop(2) = uicontrol('Style','Popup','pos',[235 110 140 20],'backgroundcolor','w',...
        'String',['9   (512 steps)|8   (256 steps)|7   (128 steps)|6   (64 steps)|5   (32 steps)',...
            '|4   (16 steps)|3   (8 steps)|2   (4 steps)|1   (2 steps)']);
    
    method(1) = uicontrol('style','text','backgroundcolor',[0.8 0.8 0.9],'pos',[45 146 170 20],...
        'horizon','right','string',['Choose a numerical method :'],'visible','off','fontweight','bold');
    
    method(2) = uicontrol('style','text','backgroundcolor',[0.8 0.8 0.9],'pos',[70 95 155 40],...
        'horizon','center','string',['Number of steps to perform before plot (in powers of 2) :'],'visible','off');
    
    fcheckLabel = uicontrol('style','text','backgroundcolor',[0.8 0.8 0.9],'pos',[45.5 70 100 20],...
        'horizon','right','string','Draw field lines :','visible','off');
    fcheckWork = ['sh = get(gcf,''user'');', 'fval = get(sh(56),''value'');',...
            'if (fval == 1) set(sh(57),''enable'',''on'',''backgroundcolor'',''w'');',...
            'else set(sh(57),''enable'',''off'',''backgroundcolor'',[0.9 0.9 0.9]); end;','set(gcf,''user'',sh);'];
    field(1) = uicontrol('style','checkbox','backgroundcolor',[0.8 0.8 0.9],'pos',[155 73 20 20],...
        'horizon','right','value',1,'visible','off','call',fcheckWork);
    feditLabel = uicontrol('style','text','backgroundcolor',[0.8 0.8 0.9],'pos',[200 70 120 20],...
        'horizon','right','string','No. of field segments :','visible','off');
    field(2) = uicontrol('style','edit','backgroundcolor','w','pos',[325 73 50 20],...
        'horizon','center','string',num2str(20),'visible','off','enable','on');
    
    quitButt =['sh = get(gcf,''user'');','tempDir = get(sh(2),''user'');','numerical(tempDir,''quit'')'];
    procButt =['sh = get(gcf,''user'');','tempDir = get(sh(2),''user'');','numerical(tempDir,''proceed'')'];
    closePlotButt =['sh = get(gcf,''user'');','tempDir = get(sh(2),''user'');','numerical(tempDir,''closeplot'')'];
    
    butt(1) = uicontrol('style','push','pos',[82 23 70 20],'fontweight','bold',...
        'string','Quit','call',quitButt,'visible','off');
    
    butt(2) = uicontrol('style','push','pos',[317 23 70 20],'fontweight','bold',...
        'string','Proceed','call',procButt,'visible','off');
    
    butt(3) = uicontrol('style','push','pos',[200 23 70 20],'fontweight','bold',...
        'string','Close Plots','call',closePlotButt,'visible','off');
    
    set(dfset,'CloseRequestFcn',quitButt);
    
    v_name = uicontrol('style','text','pos',[400 397 70 15],'foregroundcolor',[0.4 0.4 0.4],...
        'string','Version 1.7','backgroundcolor',[0.95 0.85 0.75]);
    
    % Record the handles in the User Data of the Set Up figure.
    
    sh = [frame,eq,wind,method,WINvect,ysoln,error,order,field,pop];
    set(gcf,'user',sh);
    set(get(gcf,'child'),'vis','on');
    
    % Record the important information in the User Data of
    % varous controls.
    sm = str2mat(Tname,derivstr);
    set(sh(1),'user',sm);
    set(sh(2),'user',workingDir);
    
    % Proceed and print sequences
elseif strcmpi(action,'proceed') | strcmpi(action,'print')
    
    % Proceed connects Setup with the Display window.
    sh       = get(gcf,'user');
    sm       = get(sh(1),'user');
    
    Tname    = deblank(sm(1,:));
    derivstr = deblank(sm(2,:));
    
    if strcmpi(action,'proceed')
        sh(25)   = get(sh(58),'Value');
        % If an old function m-file exists delete it, and then build a new one.
        tee = clock; tee = ceil(tee(6)*100);
        FFFCCNNN=['dftp',num2str(tee)];
        
        fcnstr = ['function YyYypr = ',FFFCCNNN,'(TtTt,YyYy)\n\n'];
    
        varstr = [Tname,' = TtTt; x = YyYy;\n\n'];
        lenstr = ['l = length(YyYy);\n'];
        derstr1 = ['YyYypr = ', derivstr,';\n'];
        derstr2 = ['if (length(YyYypr) < l) YyYypr = YyYypr*ones(1,l);end\n'];
        
        dff = fopen([FFFCCNNN,'.m'],'w');
        fprintf(dff,fcnstr);
        fprintf(dff,varstr);
        fprintf(dff,lenstr);
        fprintf(dff,derstr1);
        fprintf(dff,derstr2);
        fclose(dff);

        fieldLine(1) = get(sh(56),'Value');
        fieldLine(2) = str2num(get(sh(57),'String'));
        noSteps  = 10 - get(sh(59),'Value');
        try
            [ysoln,error,order]=numsolve(FFFCCNNN,sh(21:25),derivstr,noSteps,Tname,fieldLine);
        catch
            delete([FFFCCNNN,'.m']);
        end
        delete([FFFCCNNN,'.m']);
        
        sh(26:35)  = ysoln;
        sh(36:45)  = error;
        sh(46:55) = order;
    end;
    
    %    Set variables are as follows -
    %    WINvect  = sh(21:25);
    %    ysoln    = sh(26:35);
    %    error    = sh(36:45);
    %    order    = sh(46:55);
    
    isWin = strncmpi(getenv('OS'),'Windows',4);
    
    % Numerical Display
    % Set up the bulletin window.
    if strcmpi(action,'proceed')
        dfdisp = figure('name','NUMERICAL Display','numb','off','visible','off','menu','none');
        figure(dfdisp);
        nframe = uicontrol('style','frame','units','normal','backgroundcolor','w','pos',[.1 .17 .8 .8],'visible','off');
    else
        dfdisp = figure('name','NUMERICAL Print Output','numb','off','visible','off','color','w','menu','none');
        figure(dfdisp);
        nframe = 0;
    end
    set(dfdisp,'CloseRequestFcn',['warndlg(cellstr({''You can only close this figure by clicking'',',...
            '''''''Close'''' button from this window or from'',',...
            '''main interface by clicking ''''Close Plots'''' button''}))']);
    if length(Tname)>1 h1str = ['Differential equation: dx/d(' Tname ') = ',derivstr];
    else h1str = ['Differential equation: dx/d' Tname ' = ',derivstr]; end;
    heading(1) = uicontrol('style','text','units','normal','pos',[.15 .84 .5 .1],'horiz','left','fontsize',12,...
        'backgroundcolor','w','string',h1str,'visible','off');
    if sh(25)==1,
        heading(2) = uicontrol('style','text','units','normal','pos',[.15 .78 .5 .1],'horiz','left',...
            'backgroundcolor','w','string','Method: Euler','fontweight','bold','visible','off');
    elseif sh(25)==2,
        heading(2) = uicontrol('style','text','units','normal','pos',[.15 .78 .5 .1],'horiz','left',...
            'backgroundcolor','w','string','Method: Improved Euler','fontweight','bold','visible','off');
    else
        heading(2) = uicontrol('style','text','units','normal','pos',[.15 .78 .5 .1],'horiz','left',...
            'backgroundcolor','w','string','Method: 4th order Runge-Kutta','fontweight','bold','visible','off');
    end;
    
    solnXcond = num2str(sh(24)); if sh(24) == inf solnXcond = 'unknown'; end
    conditions = uicontrol('style','text','units','normal','pos',[.15 .73 .74 .1],'horiz','left','backgroundcolor','w',...
        'string',['Initial ' Tname ' = ' num2str(sh(21)),'     Final ', Tname, ' = ' num2str(sh(22)),...
            '     Initial x = ',num2str(sh(23)),'     Solution at final ',Tname,' = ',solnXcond],'visible','off');
    
    if strcmpi(action,'proceed')
        % Set up the buttons.
        closeButt =['sh = get(gcf,''user'');','tempDir = get(sh(2),''user'');','numerical(tempDir,''close'',gcf)'];
        printButt =['sh = get(gcf,''user'');','tempDir = get(sh(2),''user'');','numerical(tempDir,''print'',gcf)'];
        dbutt(1) = uicontrol('style','push','units','normal','pos',...
            [0.6 0.07 0.1 0.06],'fontweight','bold','string','Close',...
            'call',closeButt,'visible','off');
        
        dbutt(2) = uicontrol('style','push','units','normal','pos',...
            [0.3 0.07 0.1 0.06],'fontweight','bold','string','Print',...
            'call',printButt,'visible','off');
    else dbutt = [0 0];
    end
    
    if sh(24)~=inf,
        for i=1:10,
            if (isWin)
                solnpos = [.27 .7-i*0.05 .2 .04];
                steppos = [.165 .7-i*0.05 .05 .04];
                errpos = [.51 .7-i*0.05 .2 .04];
                orrpos = [.75 .65-i*0.05 .13 .04];
            else
                solnpos = [.185 .7-i*0.05 .3 .04];
                steppos = [.185 .7-i*0.05 .05 .04];
                errpos = [.48 .7-i*0.05 .25 .04];
                orrpos = [.76 .65-i*0.05 .13 .04];
            end
            soln(i)=uicontrol('style','text','units','normal','pos',solnpos,...
                'horiz','right','backgroundcolor','w','string',num2str(sh(25+i),'%0.12e'),'visible','off');  
            step(i)=uicontrol('style','text','units','normal','pos',steppos,...
                'horiz','right','backgroundcolor','w','string',num2str(2^(i-1)),'visible','off');  
            err(i)=uicontrol('style','text','units','normal','pos',errpos,...
                'backgroundcolor','w','horiz','right','string',num2str(sh(35+i),'%0.12e'),'visible','off');
            if i<10
                orstr = num2str(sh(45+i),'%0.10f');
                for k = 1:4
                    if sh(45+i) == k*j orstr = ['NaN(' num2str(k) ')']; break; end;
                end
                orr(i)=uicontrol('style','text','units','normal','pos',orrpos,...
                    'horiz','right','backgroundcolor','w','string',orstr,'visible','off');
            end
        end;
        
        if (isWin)
            head3pos = [.11 .72 .3 .04];
            head4pos = [.33 .72 .3 .04];
            head5pos = [.672 .72 .1 .04];
            head6pos = [.76 .72 .135 .04];
        else
            head3pos = [.105 .72 .3 .04];
            head4pos = [.32 .72 .3 .04];
            head5pos = [.67 .72 .1 .04];
            head6pos = [.745 .72 .15 .04];
        end
        
        heading(3) = uicontrol('style','text','units','normal','pos',head3pos,...
            'horiz','left','backgroundcolor','w','string','No. of Steps','visible','off');
        heading(4) = uicontrol('style','text','units','normal','pos',head4pos,...
            'horiz','left','backgroundcolor','w','string',['Solution at final ' Tname],'visible','off'); 
        heading(5) = uicontrol('style','text','units','normal','pos',head5pos,...
            'horiz','left','backgroundcolor','w','string','Error','visible','off');
        heading(6) = uicontrol('style','text','units','normal','pos',head6pos,...
            'horiz','left','backgroundcolor','w','string','Effective order','visible','off');
        if strcmpi(action,'proceed') for i = 3:6 set(heading(i),'foregroundcolor','b'); end; end;
        set([nframe,err,orr,heading,dbutt,step,soln,conditions],'visible','on');
        
    else
        
        for i=1:10,
            if (isWin)
                solnpos = [.35 .7-i*0.05 .2 .04];
                steppos = [.15 .7-i*0.05 .1 .04];
                orrpos = [.58 .6-i*0.05 .2 .04];
            else
                solnpos = [.265 .7-i*0.05 .3 .04];
                steppos = [.175 .7-i*0.05 .1 .04];
                orrpos = [.605 .6-i*0.05 .2 .04];
            end
            soln(i)=uicontrol('style','text','units','normal','pos',solnpos,...
                'horiz','right','backgroundcolor','w','string',num2str(sh(25+i),'%0.12e'),'visible','off');  
            step(i)=uicontrol('style','text','units','normal','pos',steppos,...
                'horiz','right','backgroundcolor','w','string',num2str(2^(i-1)),'visible','off');  
            if i<9
                orstr = num2str(sh(45+i),'%0.10f');
                for k = 1:4
                    if sh(45+i) == k*j orstr = ['NaN(' num2str(k) ')']; break; end;
                end
                orr(i)=uicontrol('style','text','units','normal','pos',orrpos,...
                    'horiz','right','backgroundcolor','w','string',orstr,'visible','off');
            end;
        end;
        
        heading(3) = uicontrol('style','text','units','normal','pos',[.15 .72 .3 .04],...
            'horiz','left','backgroundcolor','w','string','No. of Steps','visible','off');
        heading(4) = uicontrol('style','text','units','normal','pos',[.38 .72 .2 .04],...
            'horiz','center','backgroundcolor','w','string',['Solution at final ' Tname],'visible','off'); 
        heading(5) = uicontrol('style','text','units','normal','pos',[.66 .72 .2 .04],...
            'horiz','left','backgroundcolor','w','string','Effective order','visible','off');
        if strcmpi(action,'proceed') for i = 3:5 set(heading(i),'foregroundcolor','b'); end; end;
        set([nframe,orr,heading,dbutt,step,soln,conditions],'visible','on');
    end;
    set(gcf,'user',sh);
    
    if strcmpi(action,'print')
        figs    = get(0,'children');
        current = length(figs);
        capture(current); 
        print(current+1);
        delete(current+1);
        delete(current);
    end
    
    % Closeplot and quit sequences
elseif strcmpi(action,'closeplot') | strcmpi(action,'quit')
    
    if strcmpi(action,'quit')
        selection = questdlg(strvcat('This will quit Numerical and close all the',...
            'corresponding figures, do you wish to continue?'),...
            'Quit Numerical','Yes','No','Yes');
    else selection = 'Yes'; end
    if strcmp(selection, 'Yes')
        figs = get(0,'children');
        for (figno = 1:length(figs))
            nn = get(figs(figno),'name');
            if strcmpi(action,'closeplot')
                TTT = findstr(nn,'NUMERICAL PLOT');
                SSS = findstr(nn,'NUMERICAL Display');
                RRR = findstr(nn,'NUMERICAL Print');
                QQQ = findstr(nn,'WARNING!!');
                if (~isempty(TTT) | ~isempty(SSS) | ~isempty(RRR) | ~isempty(QQQ)) delete(figs(figno)); end;
            elseif strcmpi(action,'quit')
                TTT = findstr(nn,'NUMERICAL');
                SSS = findstr(nn,'Error');
                RRR = findstr(nn,'WARNING!!');
                if (~isempty(TTT) | ~isempty(SSS) | ~isempty(RRR)) delete(figs(figno)); end
            end
        end
    end
    
    % Close sequence for the closure of particular output
elseif strcmpi(action,'close')
    
    figno = input1;
    delete(input1);
    if strcmp(get(figno-1,'name'),'NUMERICAL PLOT')
        delete(input1-1);
    end
    
end;