Below you will find some Matlab scripts I wrote many years ago (2001-2002), for a postgraduate course on Numerical Analysis. These scripts solve twenty or so trivial problems using (and demonstrating) various numerical methods such as: Monte Carlo integration, uniform random deviates, Runge-Kutta with Bisection, Runge Kutta with Finite-Differences Newton-Raphson, Wilkinson deflation, Broyden-Charles and others.

Many of these scripts, are obsolete or deprecated, but it’s very easy to adapt them. Most of them can also be executed in Octave, with minor alterations…

 
      % This Script finds the machine's precision using plain arithmetic.
      %
      % <- Dimitris B. Kalamaras  ->
 
      clear all;
      clc;
 
      disp( 'Machine''s precision');
      tol=1.00;
      e=1.00+tol;
      while ( e > 1.0000000 )
         tol=tol/2.00;
         e=1.00+tol;
      end
      disp(' ');
      disp('TOL of the machine is: ');
      disp(tol);
 
 
 
 
 
 
      %
      % Script for the approximation via Continuous Fraction Expansion of 2^sqrt(2)
      % 
      % It seems that the continued fraction expansion is smaller and then bigger
      % and then smaller again and then bigger again and so on....
      %
      % <- Kalamaras Dimitris ->
 
      clc;
      echo off;
      clear all;
 
      disp( 'Approximation via Continuous Fraction Expansion of 2^sqrt(2)');
 
      n=input ('Enter Expansion Depth (1..40): ');
 
      format long;
      r=zeros(n,1);
      q=zeros(n,1);
      a=zeros(n,1);
 
      a(1)=vpa(2^sqrt(2));         % Enter another number if you wish...to approximate
 
      for k=1:n 
         q(k)=fix ( a(k) );   %finds the integer part of a(k)
         r(k)=  a(k)-q(k)  ;    %finds the remainder
         a(k+1)= 1/ r(k);     %inverts the remainder, which is a(k+1)
      end
 
      for k=1:n
         number = help_function (1,k,q,r);
        if (double(a(1))-double(number))> 0 
           sprintf('Number No %d is %30.15f',k, number)
           sprintf('and is SMALLER than 2^sqrt(2) = %30.15f', a(1))
        elseif (double(a(1))-double(number))< 0
           sprintf('Number No %d is %30.15f',k,number)
           sprintf('and is BIGGER than 2^sqrt(2) = %30.15f', a(1))
        else
         sprintf('This machine reached its precision limits: 2^sqrt(2) = %30.15f', number )
         end
      end
      format;
 
 
 
 
            function y = help_function (i,n,q,r)	
            %
            %
            % The "help_function" is used for the computation of the Continued Fraction
            % 
            %
            % This "help_function" uses the idea of recursion. 
            % A function that calls itself is said to be recursive. 
            %
            % <- Dimitris Kalamaras ->
 
            % clc;
 
            if nargin==0 error('Usage: help_function(i,n,q,r), where n integer'); end
 
            if i < n
                y=q(i)+1/help_function(i+1,n,q,r);
               else 
                  y=q(n);
            end
 
 
 
 
 
 
 
      %
      % This script uses the Newton Raphson method in order to find the solutions
      % of the equation :
      % cos(x)-exp(-x^2)=0
      %
      % The scripts demands an estimation x0 and returns the approximate root.
      %
      %
      % <-Kalamaras Dimitris->
 
      clc;
      clear all;
      format long;
      %x=zeros(140,1);
 
      disp(' NEWTON-RAPHSON method the equation : ');
      disp('  cos(x)-exp(-x^2)=0 ');
 
      x(1)=input('Please enter an estimation of the root x0 : ');
      repeats=input('Enter amount of repeats for the Newton-Raphson Method : ');
      tol=input('Enter desired precision :');
 
      for k=1:repeats
         f=[cos(x(k))-exp( -x(k)*x(k) )];
         df=[-sin(x(k))+2*x(k)*exp( -x(k)*x(k) )];
 
         x(k+1)=x(k)-f/df ;
 
        if abs( x(k+1)-x(k) )< tol 
           root=x(k+1);
      	  sprintf('The root is : %20.10f found with %d repeats', root, k+1)
           break
        end
      end  
 
      if abs( x(k+1)-x(k) )> tol 
         sprintf('Could not converge to a solution with the given precision...!')
      end
 
 
 
 
 
 
 
 
      % This "nr" function uses Newton Raphson method for finding the roots of a
      % given function f
      %
      % <- Kalamaras Dimitris ->
      %
      clear all;
      clc;
      format long;
 
      x=zeros(40,1);
 
      disp('Newton Raphson method for the equation cos(x)-exp(-0.5x)');
 
      x(1)=input('Enter guess value x(1) :');
 
      %f=[1 0 0 -1] 
      %df=polyder(f)
      %r=roots(f)
 
      for k=1:40
         f=cos(x(k))-exp(-0.5*x(k));
      df=-sin(x(k))-0.5*exp(-0.5*x(k));
 
       x(k+1)=x(k)+f/df;
        if (x(k+1)-x(k))< 0.001 
           break 
        end
      end  
 
      root=x(k+1);  
      disp('The solution of cosx-exp(-0.5x)=0 in the given interval is :');
      disp(root);   
 
 
 
 
 
 
      %
      % This script uses the Newton Raphson method in order to find the solutions
      % of the equation :
      % cos(x)-exp(-x^2)=0
      %
      % The scripts demands an estimation x0 and returns the approximate root.
      %
      %
      % <-Kalamaras Dimitris->
 
      clc;
      clear all;
      format long;
      %x=zeros(140,1);
 
      disp(' NEWTON-RAPHSON with FINITE DIFFERENCES method for the equation : ');
      disp('  cos(x)-exp(-x^2)=0 ');
 
      x(1)=input('Please enter an estimation of the root x0 : ');
      repeats=input('Enter amount of repeats for the Newton-Raphson Method : ');
      tol=input('Enter desired precision :');
      h=0.005; 
      dummy=input('Starting with h = 0.005. PRESS  ');
      for i=1:10
      for k=1:repeats
         f=[cos(x(k))-exp( -x(k)*x(k) )];
 
         fh= [cos(x(k)-h)-exp( -(x(k)-h)^2 )];
 
      	df=(f-fh)/h;
 
         x(k+1)=x(k)-f/df ;
 
        if abs( x(k+1)-x(k) )< tol 
           root=x(k+1);
      	  sprintf('The root is : %20.10f found with %d repeats', root, k+1)
           break
        end
      end  
 
      if abs( x(k+1)-x(k) )> tol 
         sprintf('Could not converge to a solution with the given precision...!')
      end
 
      h=h*2;
      sprintf('Now step is h= %6.5f. PRESS  ', h)
      dummy=input (' ' );
      end
 
 
 
 
 
      % This script is an implemention of the Bisection (Bolzano)
      % method which finds a root of a function, manifestated in f_bis.m,
      % in a given space (a,b) (explicitly).
      %
      % <- Kalamaras Dimitris ->
      %
      clear all;
      clc;
 
      format long;
 
      disp('BISECTION method for the equation cos(x)-exp(-x^2)=0');
      a=input('Enter left interval bound a : ');
      b=input('Enter right interval bound b : ');
      repeats=input('Enter amount of repeats : ');
      tol=input('Enter desired precision : ');
      x=zeros(41,1);
      x0=a;
 
      x(1)=x0+((b-a)/2^1)*sign(f_bis(a))*sign(f_bis(x0));
 
      for k=1:repeats
         x(k+1)=x(k)+((b-a)/2^(k+1))*sign(f_bis(a))*sign(f_bis(x(k)));
 
         if abs(x(k+1)-x(k))< tol 
           root=x(k+1);	
      	  sprintf('The root in the given interval is : %20.15f in %d times', root, k+1) 
           break ;
        end
 
      end
 
      if abs( x(k+1)-x(k) ) > tol 
         sprintf('Could not converge to a solution with the given precision ')
      end
 
 
 
            function y=f_bis(x)
            %  function 'f_bisn' is the the function 
            %  used in the 'bisect' method
            %
            % <- Kalamaras Dimitris ->
 
            global c;
 
            if nargin==0 x=0:.1:10; end
            %if nargin==0 error('Usage: f_bis( x ) , where x whatever.'); end
 
            y=cos(x)-exp(-x*x);
 
            % with some matlab native noise!
 
            %noise=-ran0*c;
            %noise=randn*c;
            %y=y+noise*y;
 
 
 
 
 
 
      %
      % This script uses the Newton Raphson method in order to find the solutions
      % of the equation :
      % cos(x)-exp(-x^2)=0
      %
      % The scripts demands an estimation x0 and returns the approximate root.
      %
      %
      % <-Kalamaras Dimitris->
      clc;
 
      x=zeros(140,1);
 
      disp(' SECANT method for the equation :');
      disp('  cos(x)-exp(-x^2)=0 ');
 
      x(1)=input('Please enter a first value x0 : ');
      x(2)=input('Please enter a first value x1 : ');
 
      repeats=input('Enter amount of repeats for the Secant Method : ');
      tol=input('Enter desired precision : ');
 
      for k=2:repeats
 
         fkp=[cos(x(k-1))-exp( -x(k-1)*x(k-1) )];
         fk=[cos(x(k))-exp( -x(k)*x(k) )];
         x(k+1)=x(k)- fk * ( x(k)-x(k-1))/(fk-fkp);
 
        if abs(x(k+1)-x(k))< tol  
           root=x(k+1);
           sprintf('The root is : %15.5f found with %d repeats', root, k+1)
           break; 
        end
      end  
 
      if abs(x(k+1)-x(k))> tol  break; 
         sprintf('Could not converge to a solution with the given precision')
      end
 
 
 
 
 
 
 
      % This script examines the behaviour of Newton Raphson  and BISECTION methods
      % concerning convergence. It applies white GAUSSIAN noise with various values
      % of variance and mean value zero for the solution of the equation :
      % cos(x)-exp(-x^2)=0
      %
      % <- Dimitris B. Kalamaras ->
 
      clc;
      clear all;
      format long;
      disp(' Convergence of NEWTON-RAPHSON & BISECTION methods for the equation : ');
      disp('  cos(x)-exp(-x^2)=0 ');
      disp(' by applying white Gaussian Noise for various sigmas');
      disp(' ');
 
      disp('NEWTON RAPHSON: estimation x(1)=5');
      x(1)=5;
      disp('NEWTON RAPHSON & BISECTION: repeats 300 max');
      repeats=300; 
      disp('NEWTON RAPHSON & BISECTION: precision 0.0001');
      tol=0.0001;
      disp('BISECTION: left interval bound a = 2 ');
      a=2;
      disp('BISECTION: right interval bound b =5 ');
      b=5;
      s=0.00001;
      key=0;
      disp(' ');
      dummy=input('The first results will have no Gaussian Noise.    PRESS ');
 
      for i=1:10
      % here starts the NEWTON RAPHSON algorithm...
      for k=1:repeats
         f=[cos(x(k))-exp( -x(k)*x(k) )]+key*gaussian(0,s);
         df=[-sin(x(k))+2*x(k)*exp( -x(k)*x(k) )];
 
         x(k+1)=x(k)-f/df ;
 
        if abs( x(k+1)-x(k) )< tol 
           root=x(k+1);
      	  sprintf('NEWTON - RAPHSON : The root is  %20.10f found with %d repeats', root, k+1)
           break 
        end
      end  
 
      if abs( x(k+1)-x(k) )> tol 
         sprintf('Could not converge to a solution with the given precision...!')
      end
 
 
      % Here starts the BISECTION algorithm...
 
      x0=a;
 
      fa=[cos(a)-exp( -a*a )]+key*gaussian(0,s);
      fx0=[cos(x0)-exp( -x0*x0 )]+key*gaussian(0,s);
 
      x(1)=x0+((b-a)/2^1)*sign(fa)*sign(fx0);
 
      for k=1:repeats
         fa=[cos(a)-exp( -a*a )]+key*gaussian(0,s);
      	fxk=[cos(x(k))-exp( -x(k)*x(k) )]+key*gaussian(0,s);
         x(k+1)=x(k)+( (b-a)/2^(k+1))*sign(fa)*sign(fxk) ;
 
         if abs(x(k+1)-x(k))< tol 
           root=x(k+1);	
      	  sprintf('BISECTION : The root in the given interval is  %20.10f in %d repeats', root, k+1) 
           break ;
        end
 
      end
 
      s=s*3;
 
      sprintf('Now Using Gaussian Noise with sigma s= %8.5f .     PRESS ',s)
      dummy=input('');
 
      clc;
      x(1)=5;
      key=1;
      end
 
 
 
      disp('Apparently, the BISECTION method endures noise for all sigma''s');
 
 
            function  Y=gaussian ( mean, sigma)
            % 
            % Function to return random variable distributed
            % according to Gaussian distribution with mean mean 
            % and standard deviation sigma.
            % It uses the "ACCEPT-REJECT" method
            %
            % < Dimitris B. Kalamaras > 
 
            if nargin < 2 error('Usage : gaussian (mean, sigma) '); end;
 
             ymin = mean - 4. * sigma; 
             ymax = mean + 4. * sigma; 
             Pymax = 1. / sqrt (2. * pi) / sigma;  % this is pymax ....
            % Calculate random value uniformly distributed 
            % in range ymin to ymax
             y = ymin + (ymax - ymin) *  rand; %/ double (RANDMAX); 
            % Calculate Py
             Py = exp (- (y - mean) * (y - mean) / 2. / sigma / sigma) / sqrt (2. * pi) / sigma;
            % Calculate random value uniformly distributed in range 0 to Pymax
             x = Pymax * rand ; %/ double (RANDMAX); 
            % If x > Py reject value and recalculate
            if x > Py 
             %	 disp(x), disp(Py), disp(y); 
              Y=gaussian (mean, sigma); 
             else  Y=y; 
            end
 
 
 
 
 
 
      %
      % This script implements the WILKINSON deflation technique using Newton Raphson method to find 
      % the solutions of the equation :
      % (x-1)^2*(x+2)*x=0
      %
      % The scripts demands assumes x0=4 and returns the approximate roots over the domain of the function.
      %
      %
      % <-Kalamaras Dimitris->
      %
 
      clc;
      clear all;
      format long;
 
      disp(' Implementation of the WILKINSON technique using NEWTON-RAPHSON method for the equation : ');
      disp('  (x-1)^2(x+2)x=0 ');
      disp('Estimation of the root:  x0 = 4 ');
 
      x(1)=4;
      repeats=input('Enter amount of repeats for the Newton-Raphson Method  ( > 50 ): ');
      tol=input('Enter desired precision [ < 0.0000001 ]:');
 
      %first loop - first root    (x0=1 multiplicity 2)
      root_number=1;
      for k=1:repeats
         f=[(x(k)-1)^2 * (x(k)+2)*x(k) ];
         df=[2*( x(k)-1 ) * (x(k)+2)* x(k) + (x(k)-1)^2 * x(k) + (x(k)-1)^2 * (x(k)+2)];
 
         x(k+1)=x(k)-f/df ;
 
        if abs( x(k+1)-x(k) )< tol 
           root(root_number)=x(k+1);
           sprintf('The %1dst root is : %20.10f found with %d repeats',root_number, root(root_number), k+1)
           root_number=root_number+1;
           for i=1:k+1      % vector of repeats used in plot "against" command....
      	   y(i)=i;
      	  end
           plot(y,x,'rx:');
            clear x; x(1)=4; 
              %second loop  - second root (x0=1 multiplicity 2)
          	 for k=1:repeats
        			 f=[ (x(k)-1) * (x(k)+2)*x(k) ];
      	  		 df=[ (x(k)+2)* x(k) + (x(k)-1) * x(k) + (x(k)-1) * (x(k)+2)];
 
      	  		 x(k+1)=x(k)-f/df ;
 
       			 if abs( x(k+1)-x(k) )< tol 
         		  root(root_number)=x(k+1);
         		  sprintf('The %1dst root is : %20.10f found with %d repeats',root_number, root(root_number), k+1)
                 root_number=root_number+1;
                 clear y; 
                   for i=1:k+1      % vector of repeats used in plot "against" command....
       	   		   y(i)=i;
         			  end
                    hold on ; plot(y,x,'bx:');
                    clear x; x(1)=4; 
                  %third loop  - third root (x0=0)
        				for k=1:repeats
      		  			 f=[ (x(k)+2)*x(k) ];
      	  				 df=[ x(k) + (x(k)+2)];
 
      	  				 x(k+1)=x(k)-f/df ;
 
       					 if abs( x(k+1)-x(k) )< tol 
         					  root(root_number)=x(k+1);
      			   		  sprintf('The %1dst root is : %20.10f found with %d repeats',root_number, root(root_number), k+1)
      				    	  root_number=root_number+1;
                             clear y;
         	                 for i=1:k+1      % vector of repeats used in plot "against" command....
      							   y(i)=i;
      							  end
                             hold on;plot(y,x,'go:');
                             clear x; x(1)=4; 
      		      	    %fourth loop  - fourth root (x0=-2)
        						for k=1:repeats
      				  			 f=[ (x(k)+2) ];
      	  						 df=[1];
 
      	  						 x(k+1)=x(k)-f/df ;
 
       							 if abs( x(k+1)-x(k) )< tol 
         						  root(root_number)=x(k+1);
      				   		  sprintf('The %1dst root is : %20.10f found with %d repeats',root_number, root(root_number), k+1)
                             root_number=root_number+1; clear y;
                             for i=1:k+1      % vector of repeats used in plot "against" command....
      							   y(i)=i;
      							  end
                             hold on;plot(y,x,'yd:'); text (10,3,'Red no deflation, Blue 1st,Green 2ond and Yellow 4th deflation.');
        	    					  break
      		 					 end
      	   			 	end  
 
      						if abs( x(k+1)-x(k) )> tol 
      					  	 sprintf('Could not converge to a solution with the given precision...!')
      						end
 
 
       					  break
      		 			 end
      	   		 end  
 
      				if abs( x(k+1)-x(k) )> tol 
      			  	 sprintf('Could not converge to a solution with the given precision...!')
      				end
 
 
       			  break
       			 end
      	    end  
 
      		if abs( x(k+1)-x(k) )> tol 
      	  	 sprintf('Could not converge to a solution with the given precision...!')
      		end
 
 
           break
        end
      end  
 
      if abs( x(k+1)-x(k) )> tol 
         sprintf('Could not converge to a solution with the given precision...!')
      end
 
 
 
 
 
 
      %
      % CONTOUR plots for the the system 
      % f1= x.^2-y.^2-1 ;
      % f2= x.^2+y.^2-2;
      % 
      %
      % <- Kalamaras Dimitris ->
 
      function contr1(x1,x2);
 
      if nargin==0 
         dx=1/8;
         [x y]=meshgrid(-3:dx:3);
       end;
 
 
 
      disp( 'CONTOUR plots for the the system  : ');
      disp( ' f1= x.^2-y.^2-1 ');
      disp( ' f2= x.^2+y.^2-2 ');
 
 
 
      f1= x.^2-y.^2-1 ;
      f2= x.^2+y.^2-2;
      mesh(x,y,f1);
      hold on;
      mesh(x,y,f2);
      title('Mesh plot for f1=x^2-y^2-1 and f2=x^2+y^2-2');
 
      figure
      surf(x,y,f1);
      hold on;
      surf(x,y,f2);
      title('Surf plot for f1=x^2-y^2-1 and f2=x^2+y^2-2');
 
      figure
      contour(x,y,f1);
      title('Contour plot for f1=x^2-y^2-1');
 
      figure
      contour(x,y,f2);
      title('Contour plot for f2=x^2+y^2-2');
 
      figure
      contour(x,y,f1);
      hold on;
      contour(x,y,f2);
      title('Contour plot for f1=x^2-y^2-1 and f2=x^2+y^2-2');
 
 
 
 
 
 
 
      %
      %  This Matlab script approximates the solution of the system 
      %  {f1==x1^2 + x2^2 - 9 = 0
      %  {f2==x1   + x2   - 3 = 0   , where f1,f2 are R^2-->R, 
      %  using the Newton Method in 2 dimension:
      %  Jf(xp)*Sp=-F(xp)
      %  x(p+1)=x(p)+Sp
      %  
      %
      % <-Kalamaras Dimitris->
      %
 
      clear all;
      clc;
 
 
      disp('NEWTON 2d for the solution of the following system:');
      disp('  f1(x1,x2)=x1^2 + x2^2 - 9');
      disp('  f2(x1,x2)=x1   + x2   - 3');
      disp('  with initial values (x1,x2)=(1,5)');
      %  f1(x1,x2)=x1^2-x2^2-9
      %  f2(x1,x2)=x1+x2-3
 
 
      x1=1;             %initial values
      x2=5;
 
      max_repeats=input('Enter maximum amount of repeats for Newton 2d :');
 
      %If you want another system, change the two next lines.
      fxr=[x1^2+x2^2-9 ; x1+x2-3];    %Functions' f1 & f2 values vector
      Jf=[2*x1 , 2*x2; 1 , 1 ];     % Jacobian of f1, f2 ref to x1,x2
 
      for k=1:max_repeats
 
      	Sp=Jf\(-fxr);    % Left division of Jfxr into -fxr uses LU factorization
      	x1p=x1;
      	x2p=x2;
      	x1=x1+Sp(1);
      	x2=x2+Sp(2);
      	if abs(x1-x1p)<0.0001 & abs(x2-x2p)<0.0001 break; end;
 
         fxr=[x1^2+x2^2-9 ; x1+x2-3];   %Functions' values NEW vector
         Jf=[2*x1 , 2*x2; 1 , 1 ];  % NEW Jacobian of f1, f2 ref to x,y
      end;
 
      disp('Solution is:');
      root=[x1;x2];
      disp(root);
 
      disp('iterations done : '); disp(k);	
 
      %toc
 
 
 
 
 
 
      %  This Matlab script approximates the solution of the system 
      %  {f1==x1^2 + x2^2 - 9 = 0
      %  {f2==x1   + x2   - 3 = 0   , where f1,f2 are R^2-->R, 
      %  
      %  using the Broyden Charles Method which is a quasi-Newton in 2 dimensions:
      %
      % A(p)*Sp=-F(x(p))
      % x(p+1)=x(p)+Sp
      % y(p)=F(x(p+1))-F(x(p))
      % A(p+1)=A(p)+((y(p)-A(p)*Sp)*Sp')/Sp'*Sp
      % 
      % Broyden-Charles approximates the Jacobian, thus avoiding the calculation, if it is possible,
      % of the derivatives. 
      %
      % Use as initial matrix A, a non singular, ie A=I. A should be converging to Jacobian
      % at the end of the Broyden-Charles. 
      % 
      % <-Dimitris B. Kalamaras ->
      %
 
 
      clear all;
      clc;
 
      disp('	BROYDEN method for the following system:');
      disp('	f1(x1,x2)=x1^2 + x2^2 - 9 ');
      disp('	f2(x1,x2)=x1   + x2   - 3 ');
      disp('	with initial values (x1,x2)=(1,5)');
      %  f1(x1,x2)=x1^2+x2^2-9
      %  f2(x1,x2)=x1+x2-3
 
 
      x1=1;             %initial values
      x2=5;
 
      max_repeats=input('Enter maximum amount of repeats for Broyden :');
 
      A=eye(2);        %Change to whatever non singular matrix 2x2.
 
 
      %If you want another system, change the next lines.
      fxp=[x1^2+x2^2-9 ; x1+x2-3];    %Functions' f1 & f2 values vector
 
      for k=1:max_repeats
 
      	Sp=A\(-fxp);    % Left division of Jfxr into -fxr uses LU factorization
      	x1p=x1;
      	x2p=x2;
      	x1=x1+Sp(1);
      	x2=x2+Sp(2);
      	if abs(x1-x1p)<0.0001 & abs(x2-x2p)<0.0001 break; end;
 
      	fxp1=[x1^2-x2^2-9 ; x1+x2-3];                %Functions' values NEW vector
 
      	x2p=fxp1-fxp;
      	A=A+((x2p-A*Sp)*Sp')/(Sp'*Sp);
 
      	fxp=fxp1;
      end;
 
      disp('The Broyden Solution is:');
      root=[x1; x2];
      disp(root);
 
      disp('iterations done : '); disp(k);	
 
      %toc
 
 
 
 
 
      %  This Matlab script approximates the solution of the system 
      %  {f1==x1^2 + x2^2 - 9 = 0
      %  {f2==x1   + x2   - 3 = 0   , where f1,f2 are R^2-->R, 
      %  using the quasi-Newton Method with finite differences in 2 dimensions:
      %  Jf(xp)*Sp=-F(xp) where jacobian is approximated
      %  x(p+1)=x(p)+Sp
      %  
      %
      % <-Kalamaras Dimitris->
 
      clear all;
      clc;
 
     disp(' Quasi-NEWTON with finite differences');
      disp('                for the solution of the following system:');
      disp('f1(x1,x2)=x1^2-x2^2-9');
      disp('f2(x1,x2)=x1+x2-3');
      disp('with initial values (x1,x2)=(1,5)');
 
 
      %if you want to solve another system , change it in the function 'f_qn' ....
 
      x1=1;             %initial values
      x2=5;
 
      max_repeats=input('Enter maximum amount of repeats for Newton 2d f.diff. :');
      h=input('Enter step h : ');
 
      theta1_f1=(f_qn(1,x1+h,x2)-f_qn(1,x1,x2))/h;
      theta2_f1=(f_qn(1,x1,x2+h)-f_qn(1,x1,x2))/h;
      theta1_f2=(f_qn(2,x1+h,x2)-f_qn(2,x1,x2))/h;
      theta2_f2=(f_qn(2,x1,x2+h)-f_qn(2,x1,x2))/h;
 
 
      fxr=[f_qn(1,x1,x2) ;f_qn(2,x1,x2) ];    %Functions' f1 & f2 values vector
      Jf=[theta1_f1 , theta2_f1; theta1_f2 , theta2_f2 ];     % Jacobian of f1, f2 ref to x1,x2
 
      for k=1:max_repeats
 
      	Sp=Jf\(-fxr);    % Left division of Jfxr into -fxr uses LU factorization
      	x1p=x1;
      	x2p=x2;
      	x1=x1+Sp(1);
      	x2=x2+Sp(2);
      	if abs(x1-x1p)<0.0001& abs(x2-x2p)<0.0001 break; end;
 
         theta1_f1=(f_qn(1,x1+h,x2)-f_qn(1,x1,x2))/h;
      	theta2_f1=(f_qn(1,x1,x2+h)-f_qn(1,x1,x2))/h;
      	theta1_f2=(f_qn(2,x1+h,x2)-f_qn(2,x1,x2))/h;
      	theta2_f2=(f_qn(2,x1,x2+h)-f_qn(2,x1,x2))/h;
 
      	fxr=[f_qn(1,x1,x2) ;f_qn(2,x1,x2) ];     %Functions' values NEW vector
      	Jf=[theta1_f1 , theta2_f1; theta1_f2 , theta2_f2 ];     % NEW Jacobian of f1, f2 ref to x1,x2
 
      end;
 
      disp('NEWTON finite differences found the Solution to be: ');
      root=[x1;x2];
      disp(root);
      disp('iterations done : '); disp(k);	
      %toc
 
 
 
 
      %
      % This script calls the 'newton21' (Newton for systems) function to solve the equation
      %  z^3-1-0 and reports if a given point 'goes to' which root of the equation.
      % 
      % It constructs the corresponding fractal, which shows the areas 
      % of convergence in the given space [-2,2], [-2i,2i].
      %
      %
      % <- Dimitris B. Kalamaras ->
 
 
      clc;
      clear all;
      tic;
 
 
      disp( 'Generation of fractal from the solution of z^3-1=0 ');
      disp( '              using Newton for systems.');
 
      i1=1; i2=1 ; i3=1;
      root1=[1;0];
      root2=[-1/2;-sqrt(3)/2];
      root3=[-1/2;sqrt(3)/2];
      max_iterations=input('Enter number of repeats for Newton2D :');
      for r=-2:0.025:2     %0.025  ==== 232 secs on celeron 333mhz
         for im=-2:0.025:2
            k=newton21(r,im,max_iterations);
            disp('With initial value');
            disp(r+im*i);
            disp('Newton2D converges to');
            disp(k);
            if abs(k-root1) < [0.1;0.1]
               red(i1)=r+im*i;
               i1=i1+1;
            elseif abs(k-root2) < [0.1;0.1] 
              green(i2)=r+im*i;
              i2=i2+1;
            elseif abs(k-root3)< [0.1;0.1]
             blue(i3)=r+im*i;
             i3=i3+1;
            end
         end
      end
 
      figure;
           box off;
           hold on;
 
           plot(blue,'b.');
           plot(red,'r.');
           plot(green,'g.');
           title('Fractal, areas of convergence in [-2,2][-2i,2i] for z^3-1=0 via Newton in 2 dimensions');
            xlabel('red=1, green=-1/2-sqrt(3)/2, blue=-1/2+sqrt(3)/2');
            hold off;
            zoom on;
        toc;
 
 
 
 
            function root=newton21 (a,b,max_iter)
            %  Newton2 approximates the solution of the system {f1==0,f2==0}, where f1,f2 are R^2-->R, 
            %  using the Newton Method in 2 dimension:
            %  Jf(xp)*Sp=-F(xp)
            %  x(p+1)=x(p)+Sp
 
            %
            % <-Kalamaras Dimitris->
 
 
            %clear all;
            clc;
            if nargin<3 error('Usage newton21(a,b,max_iter) where a,b rational Re(z0) & Im(z0), max_iter integer'); end;
 
            % for z^3-1=0 it is
            %  z^3-1=0
            % (x-y*i)^3-1=0
            %  x^3-3*i*x^2*y-3*x*y^2+i*y^3-1=0
            %  f1=Re(z^3-1)=x^3-3*x*y^2-1
            %  f2=Im(z^3-1)=y^3-3*x^2*y
            %tic
 
 
            x=a;             %initial values
            y=b;
 
 
            %If you want another system, change the two next lines.
 
            fxr=[x^3-3*x*y^2-1; y^3-3*x^2*y;];             %Functions' f1 & f2 values vector
 
            Jf=[3*x^2-3*y^2, -6*x*y ; -6*x*y, 3*y^2-3*x^2];     % Jacobian of f1, f2 ref to x,y
 
 
            for iteration=1:max_iter
 
            	Sp=Jf\(-fxr);    % Left division of Jfxr into -fxr uses LU factorization
            	xp=x;
            	yp=y;
            	x=x+Sp(1);
            	y=y+Sp(2);
            	if abs(x-xp)<0.0001 & abs(y-yp)<0.0001 break; end;
 
            	fxr=[x^3-3*x*y^2-1; y^3-3*x^2*y];             %Functions' values NEW vector
 
            	Jf=[3*x^2-3*y^2, -6*x*y ; -6*x*y, 3*y^2-3*x^2];     % NEW Jacobian of f1, f2 ref to x,y
 
 
            end;
 
            %disp('Solution is:');
            Jf
            iteration
            root=[x;y];
 
            %toc
 
 
 
      %
      % This script calls the 'newton21_noise' (Newton for systems with gaussian noise) function to solve the equation
      %  z^3-1-0 and reports if a given point 'goes to' which root of the equation.
      % 
      % It constructs the corresponding fractal, which shows the areas 
      % of convergence in the given space [-2,2], [-2i,2i].
      %
      %
      % <- Dimitris B. Kalamaras ->       
      %
 
      clc;
      clear all;
      tic;
 
 
      disp( 'Generation of fractal from the solution of z^3-1=0 ');
      disp( '              using Newton for systems with gaussian noise N(0,s).');
      i1=1; i2=1 ; i3=1;
      root1=[1;0];
      root2=[-1/2;-sqrt(3)/2];
      root3=[-1/2;sqrt(3)/2];
      max_iterations=input('Enter number of repeats for Newton2D :');
      s=input('Enter value of parameter sigma for the gaussian noise : ');
      for r=-2:0.025:2     %0.025  ==== 232 secs on celeron 333mhz
         for im=-2:0.025:2
            k=newton21_noise(r,im,max_iterations,s);    %calls newton21 with noise to return a point...
            disp('With initial value');
            disp(r+im*i);
            disp('Newton2D converges to');
            disp(k);
            if abs(k-root1) < [0.1;0.1]
               red(i1)=r+im*i;
               i1=i1+1;
            elseif abs(k-root2) < [0.1;0.1] 
              green(i2)=r+im*i;
              i2=i2+1;
            elseif abs(k-root3)< [0.1;0.1]
             blue(i3)=r+im*i;
             i3=i3+1;
            end
         end
      end
 
      figure;
           box off;
           hold on;
 
           plot(blue,'b.');
           plot(red,'r.');
           plot(green,'g.');
           title('Fractal- areas of convergence in [-2,2][-2i,2i] for z^3-1=0 via Newton 2d with Noise s=0.1');
            xlabel('red=1, green=-1/2-sqrt(3)/2, blue=-1/2+sqrt(3)/2');
            hold off;
            zoom on;
        toc;
 
 
 
 
 
 
            function root=newton21_noise (a,b,max_iter, s)
            %  Newton2 approximates the solution of the system {f1==0,f2==0}, where f1,f2 are R^2-->R, 
            %  using the Newton Method in 2 dimension:
            %  Jf(xp)*Sp=-F(xp)
            %  x(p+1)=x(p)+Sp
 
            %
            % <-Kalamaras Dimitris->
 
 
            %clear all;
            clc;
            if nargin<3 error('Usage newton21_noise(a,b,max_iter,s) where a,b rational Re(z0) & Im(z0), max_iter integer, s: sigma'); end;
 
            % for z^3-1=0 it is
            %  z^3-1=0
            % (x-y*i)^3-1=0
            %  x^3-3*i*x^2*y-3*x*y^2+i*y^3-1=0
            %  f1=Re(z^3-1)=x^3-3*x*y^2-1
            %  f2=Im(z^3-1)=y^3-3*x^2*y
 
 
            	x=a;             %initial values
            	y=b;
 
 
            	%If you want another system, change the two next lines.
 
            	fxr=[x^3-3*x*y^2-1+gaussian(0,s); y^3-3*x^2*y+gaussian(0,s);];       %Functions' f1 & f2 values vector
 
            	Jf=[3*x^2-3*y^2, -6*x*y ; -6*x*y, 3*y^2-3*x^2];     % Jacobian of f1, f2 ref to x,y
 
 
            	for iteration=1:max_iter
 
            		Sp=Jf\(-fxr);    % Left division of Jfxr into -fxr uses LU factorization
            		xp=x;
            		yp=y;
            		x=x+Sp(1);
            		y=y+Sp(2);
            		if abs(x-xp)<0.0001 & abs(y-yp)<0.0001 break; end;
 
            		fxr=[x^3-3*x*y^2-1; y^3-3*x^2*y];             %Functions' values NEW vector
 
            		Jf=[3*x^2-3*y^2, -6*x*y ; -6*x*y, 3*y^2-3*x^2];     % NEW Jacobian of f1, f2 ref to x,y
 
            	end;
 
            	%disp('Solution is:');
            	Jf
            	iteration
            	root=[x;y];
 
 
 
 
 
      %
      % ITERATIVE METHODS
      %
      % JACOBI - GAUSS SEIDEL - S.O.R. METHODS FOR SOLVING LINEAR NXN SYSTEMS
      %
      %  
      % <-Dimitris B. Kalamaras ->
      %
 
 
      clc;
      clear all;
 
 
      disp( 'LINEAR SYSTEMS SOLVE METHODS (iterative)');
      disp('                 JACOBI - GAUSS SEIDEL - S.O.R. methods for a NxN linear system ');
      disp('-------------------------------------------------------------------------------- ');
 
      disp( ' System''s dimension, N   ');
      n=4; disp(n);
      disp( 'Matrix A, ');
      A(1,1)=4; A(1,2)=-1; A(1,3)=-1; A(1,4)=0; 
      A(2,1)=-1; A(2,2)=4; A(2,3)=0; A(2,4)=-1;
      A(3,1)=-1; A(3,2)=0; A(3,3)=4; A(3,4)=-1;
      A(4,1)=0; A(4,2)=-1; A(4,3)=-1; A(4,4)=4;
 
      disp(A);
 
      disp( 'Matrix B,');
      B(1)=0; B(2)=0; B(3)=1000; B(4)=1000;
      disp(B);
      disp( 'Matrix X,');
      X(1)=1; X(2)=1; X(3)=1 ; X(4)=1;
      disp(X);
 
      max_repeats=15;
      disp('Maximum repeats for all three methods, ');
      disp(max_repeats);
      OMEGA=1.062;
      disp('relaxation parameter OMEGA for S.O.R. : ');
      disp(OMEGA);
      TOL=0.0001;
      disp('precision TOL : ');
      disp(TOL);
      dummie=input('Press  to Continue>');
 
      %main JACOBI-GAUSS SEIDEL- S.O.R. algorithms       (repeated for algorithmical elegance)
 
      for method=1:3
 
      	for repeats=1:max_repeats 
 
              for k=1:n
      	       previous_X(k)=X(k);
      	     end
 
      		  for i=1:n 
 
      	   	   sum1=0; sum2=0;
 
        		       if i~=1
      			   	 for j=1:i-1
      	   		   	 sum1=sum1+A(i,j)*X(j);
      		   	  	  end
          	   	 end;
 
      		       if i~=n
      	  		 		 for j=i+1:n
      			      	 sum2=sum2+A(i,j)*X(j);
         	   	     end
        	    	    end;
 
        	    	    if method==1 
      	                 Z(i)=(1/A(i,i))*(B(i)-sum1-sum2);
        	      	 elseif method==2
             	       	 X(i)=(1/A(i,i))*(B(i)-sum1-sum2);
             	    elseif method==3
              		    	 X(i)=(OMEGA/A(i,i))*(B(i)-sum1-sum2)-(OMEGA-1)*X(i);
               	 end;
              end;
 
             if method==1
             	for i=1:n
        	 		 X(i)=Z(i);
               end
             end;
 
            precision_flag=0;
             for i=1:n
                if (X(i)-previous_X(i) <= TOL) precision_flag=precision_flag+1; end
             end
             if (precision_flag==n) break; end
 
 
         end
 
        if method==1 
           disp('');
           disp('The JACOBI solution is X = ');
           X
           disp('Iterations needed : ');
           disp(repeats);
         elseif method==2
          	disp('');
         	disp('The GAUSS -SEIDEL solution is X = ');
         	X
      	   disp('Iterations needed : ');
      	   disp(repeats);
         elseif method==3
       	   disp('');
         	disp('The SOR solution is X = ');
      	   X
         	disp('Iterations needed : ');
      	   disp(repeats);
         end
 
 
      end
 
      % Under the musical influence of:
      % "The man who sold the world " (David Bowie)
 
 
 
 
 
 
 
      %
      % ITERATIVE METHODS
      %
      % JACOBI METHOD FOR SOLVING LINEAR NXN SYSTEMS
      %
      %  
      % <-Dimitris B. Kalamaras ->
 
      clc;
      clear all;
 
 
      disp( 'LINEAR SYSTEMS SOLVE METHODS (iterative)');
      disp('					  JACOBI method for a NxN linear system ');
 
 
      n=input( 'ENTER system dimension  N:');
 
      disp( 'ENTER matrix A, element by element in a row by row fashion:');
 
      for i=1:n
         for j=1:n 
            A(i,j)=input ('element : ');
      	end
      end
 
      disp( 'ENTER matrix B:');
      for  i=1:n
         B(i)=input('B element : ');
      end
 
      disp( 'ENTER matrix X:');
      for  i=1:n
         X(i)=input('X element : ');
      end
 
      max_repeats=input('ENTER maximum repeats : ');
 
 
      disp( 'MATRIX A:');
      A
      disp('MATRIX B:');
      B
 
      %main algorithm       
      for repeats=1:max_repeats 
      	 for i=1:n 
             sum1=0; sum2=0;
             if i~=1
      	   	 for j=1:i-1
         	   	 sum1=sum1+A(i,j)*X(j);
      	       end
             end
 
             if i~=n
         		 for j=i+1:n
      	      	 sum2=sum2+A(i,j)*X(j);
                end
             end
 
      	    Z(i)=(1/A(i,i))*(B(i)-sum1-sum2);
      	 end
 
       	for i=1:n
         	 X(i)=Z(i);
      	 end
 
      end
 
 
 
      disp('');
      disp('The solution is X = ');
      X
 
 
 
      %
      % ITERATIVE METHODS
      %
      % GAUSS SEIDEL METHOD FOR SOLVING LINEAR NXN SYSTEMS
      %
      %  
      % <-Dimitris B. Kalamaras ->
      %
 
 
      clc;
      clear all;
 
 
      disp( 'LINEAR SYSTEMS SOLVE METHODS (iterative)');
      disp('					  GAUSS SEIDEL method for a NxN linear system ');
 
      n=input( 'ENTER system dimension  N:');
 
      disp( 'ENTER matrix A, element by element in a row by row fashion:');
 
      for i=1:n
         for j=1:n 
            A(i,j)=input ('element : ');
      	end
      end
 
      disp( 'ENTER matrix B:');
      for  i=1:n
         B(i)=input('B element : ');
      end
 
      disp( 'ENTER matrix X:');
      for  i=1:n
         X(i)=input('X element : ');
      end
      max_repeats=input('ENTER maximum repeats : ');
 
 
      disp( 'MATRIX A:');
      A
      disp('MATRIX B:');
      B
 
      %main algorithm       
      for repeats=1:max_repeats
 
       	for i=1:n 
             sum1=0; sum2=0;
             if i ~= 1 
      	    	for j=1:i-1
      	      	 sum1=sum1+A(i,j)*X(j);
      	       end
         	 end
         	if i ~= n  
      	     for j=i+1:n
      	        sum2=sum2+A(i,j)*X(j);
              end
            end
          	 X(i)=(1/A(i,i))*(B(i)-sum1-sum2);
      	 end
 
       end
 
 
 
      disp('');
      disp('The solution is X = ');
      X
 
 
 
 
 
      %
      % ITERATIVE METHODS
      %
      % SUCCESIVE OVERRELAXATION METHOD FOR SOLVING LINEAR NXN SYSTEMS
      %
      % BEST OMEGA=1.062 :-)
      %
      % <-Dimitris B. Kalamaras ->
 
      clc;
      clear all;
 
 
      disp( 'LINEAR SYSTEMS SOLVE METHODS (iterative)');
      disp('                 Succesive Over Relaxation method for a NxN linear system ');
 
      n=input( 'ENTER system dimension  N:');
 
      disp( 'ENTER matrix A, element by element in a row by row fashion:');
 
      for i=1:n
         for j=1:n 
            A(i,j)=input ('element : ');
      	end
      end
 
      disp( 'ENTER matrix B:');
      for  i=1:n
         B(i)=input('B element : ');
      end
 
      disp( 'ENTER matrix X:');
      for  i=1:n
         X(i)=input('X element : ');
      end
      max_repeats=input('ENTER maximum repeats : ');
      OMEGA=input('ENTER relaxation parameter omega : ');
 
 
      disp( 'MATRIX A:');
      A
      disp('MATRIX B:');
      B
 
      %main algorithm       
      for repeats=1:max_repeats
 
       	for i=1:n 
             sum1=0; sum2=0;
             if i ~= 1 
      	    	for j=1:i-1
      	      	 sum1=sum1+A(i,j)*X(j);
      	       end
         	 end
         	if i ~= n  
      	     for j=i+1:n
      	        sum2=sum2+A(i,j)*X(j);
              end
            end
          	 X(i)=(OMEGA/A(i,i))*(B(i)-sum1-sum2)-(OMEGA-1)*X(i);
      	 end
 
       end
 
 
      disp('');
      disp('The solution is X = ');
      X
 
 
 
 
 
      %  steepest descent with armijo  -- Check it . I do not remember if it worked out!!!
      %
      % <- Dimitris B. Kalamaras  ->
 
      clc;
      clear all;
      format long;
      x(1,1)=input('Enter initial x1(0) [1] : ');
      x(2,1)=input('Enter initial x2(0) [1]: ');
      max_repeats=input('Enter maximum repeats [ > 1500 ]:');
      tol=input('Enter precision level tol [ < 0.00001 ]: ');
 
 
      h0=0.99;   % this has to be less than 1
 
      for m=1:45  % play with this if you like but 50 is the maximum matlab knows...heheheh
         hta(m)=h0*2^(1-m);
      end
 
      h=h0;
 
      for i=2:max_repeats
 
        % anadelta_f=[2-2*x(1,i-1) ; 2-2*x(2,i-1)];
 
         anadelta_f = [2*x(1,i-1) ; 2*x(2,i-1)];
         if anadelta_f==0 disp('zero gradient. possible minimum' ); break; end
 
         x(1,i) = x(1,i-1) - h*anadelta_f(1);
         x(2,i) = x(2,i-1) - h*anadelta_f(2);
 
         f_old = x(1,i-1)^2+x(2,i-1)^2;
         f_new = x(1,i)^2+x(2,i)^2;
         %f_old = 2+ 2*x(1,i-1)+ 2*x(2,i-1)- x(1,i-1)^2 - x(2,i-1)^2;
         %f_new = 2+ 2*x(1,i)+ 2*x(2,i)- x(1,i)^2 - x(2,i)^2;
 
 
           for m=2:10
            if (f_new-f_old) <= -hta(m)/2 * norm(anadelta_f,2)^2 h=hta(m);   end
         end
 
         if abs(x(1,i)-x(1,i-1))
 
 
 
 
 
 
      %
      % LOTKA VOLTERA 2 POPULATION MODEL
      %
      %  y'=f(x,y,z) = ay-yz
      %  z'=g(x,y,z) = -bz+yz
      %  y0=y(x0)
      %  z0=z(x0)
      %
      %  Using EULER, EULER-CAUCHY, IMPROVED-EULER, RUNGE-KUTTA 4th CLASS numerical solution methods.
      %
      %	% if you want to solve another system of ODEs change the f_ode 
      %
      % <- Dimitris B. Kalamaras  ->
 
 
      clc;
      clear all;
 
       disp('LOTKA VOLTERA 2 POPULATION MODEL');
 
      disp('');
      x(1)=input('Enter x0 [0]: '); 
      y(1)=input('Enter y0 [100]: '); 
      z(1)=input('Enter z0 [15]: '); 
      h=input ('Enter h  [0.0005]: '); 
      n=input ('Enter N  [3000]: '); 
      a=input ('Enter a  [60]: '); 
      b=input ('Enter b  [34]: '); 
 
      disp('Select a numerical solution ODE method:');
      disp('1. Euler');
      disp('2. Euler_Cauchy');
      disp('3. Improved Euler');
      disp('4. Runge 4th Class');
      choice=input('> ');
 
 
      switch (choice) 
      case 1, 
         for i=1:n
         y(i+1)=y(i)+h*f_ode(x(i),y(i),z(i), 1, a,b);
         z(i+1)=z(i)+h*f_ode(x(i),y(i),z(i), 2,a,b);
         x(i+1)=x(i)+h;
      	sprintf('x = %15.10f,  y = %15.10f,   z = %15.10f',x(i+1),y(i+1),z(i+1));
      end
      hold on;
      plot(x,y,'r*-');
      plot(x,z,'b.:');
      hold off;
      figure;
      plot(y,z,'g:');
      title('LOTKA VOLTERA 2 populations MODEL with EULER method');
 
      case 2,
 
      for i=1:n
         y(i+1)=y(i)+h*f_ode( x(i)+h/2, y(i)+(h/2)*f_ode(x(i),y(i),z(i),1,a,b), z(i)+(h/2)*f_ode(x(i),y(i),z(i),2,a,b), 1, a,b );
         z(i+1)=z(i)+h*f_ode( x(i)+h/2, y(i)+(h/2)*f_ode(x(i),y(i),z(i),1,a,b), z(i)+(h/2)*f_ode(x(i),y(i),z(i),2,a,b), 2, a,b );
         x(i+1)=x(i)+h;
         sprintf('x = %15.10f,  y = %15.10f,   z = %15.10f',x(i+1),y(i+1),z(i+1));
 
      end
      hold on;
      plot(x,y,'r*-');
      plot(x,z,'b.:');
      hold off;
      figure;
      plot(y,z,'g:');
      title('LOTKA VOLTERA 2 populations MODEL with EULER-CAUCHY method');
 
      case 3, 
         for i=1:n
       	y(i+1)=y(i)+ (h/2) * ( f_ode( x(i),y(i),z(i),1,a,b ) + f_ode( x(i)+h, y(i)+h*f_ode(x(i),y(i),z(i),1,a,b), z(i)+h*f_ode(x(i),y(i),z(i),2,a,b), 1 , a,b ) );
         z(i+1)=z(i)+ (h/2) * ( f_ode( x(i),y(i),z(i),2,a,b ) + f_ode( x(i)+h, y(i)+h*f_ode(x(i),y(i),z(i),1,a,b), z(i)+h*f_ode(x(i),y(i),z(i),2,a,b), 2 , a ,b ) );
         x(i+1)=x(i)+h;
         sprintf('x = %15.10f,  y = %15.10f,   z = %15.10f',x(i+1),y(i+1),z(i+1));
        end
      hold on;
      plot(x,y,'r*-');
      plot(x,z,'b.:');
      hold off;
      figure;
      plot(y,z,'g:');
      title('LOTKA VOLTERA 2 populations MODEL with IMPROVED EULER');
 
      case 4, 
         for i=1:n
          k1=h*f_ode(x(i),y(i),z(i),1,a,b);
          l1=h*f_ode(x(i),y(i),z(i),2,a,b);
          k2=h*f_ode(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,1,a,b);
          l2=h*f_ode(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,2,a,b);
          k3=h*f_ode(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,1,a,b);
          l3=h*f_ode(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,2,a,b);
          k4=h*f_ode(x(i)+h, y(i)+k3, z(i)+l3,1,a,b);
          l4=h*f_ode(x(i)+h, y(i)+k3, z(i)+l3,2,a,b);
 
          y(i+1)=y(i)+ (1/6)*(k1 + 2*k2 + 2*k3 + k4);
          z(i+1)=z(i)+ (1/6)*(l1 + 2*l2 + 2*l3 + l4);
 
         x(i+1)=x(i)+h;
         sprintf('x = %15.10f,  y = %15.10f,   z = %15.10f',x(i+1),y(i+1),z(i+1));    %exclamation mark here to abandon output.
      	end
      hold on;
      plot(x,y,'r*-');
      plot(x,z,'b.:');
      hold off;
      figure;
      plot(y,z,'g:');
      title('LOTKA VOLTERA 2 populations MODEL with RUNGE KUTTA 4TH CLASS');
 
      otherwise 
         for i=1:n
          k1=h*f_ode(x(i),y(i),z(i),1,a,b);
          l1=h*f_ode(x(i),y(i),z(i),2,a,b);
          k2=h*f_ode(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,1,a,b);
          l2=h*f_ode(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,2,a,b);
          k3=h*f_ode(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,1,a,b);
          l3=h*f_ode(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,2,a,b);
          k4=h*f_ode(x(i)+h, y(i)+k3, z(i)+l3,1,a,b);
          l4=h*f_ode(x(i)+h, y(i)+k3, z(i)+l3,2,a,b);
 
          y(i+1)=y(i)+ (1/6)*(k1 + 2*k2 + 2*k3 + k4);
          z(i+1)=z(i)+ (1/6)*(l1 + 2*l2 + 2*l3 + l4);
 
         x(i+1)=x(i)+h;
         sprintf('x = %15.10f,  y = %15.10f,   z = %15.10f',x(i+1),y(i+1),z(i+1));    %exclamation mark here to abandon output.
        end
      hold on;
      plot(x,y,'r*-');
      plot(x,z,'b.:');
      hold off;
      figure;
      plot(y,z,'g:');
      title('LOTKA VOLTERA 2 populations MODEL with RK 4TH CLASS ');
 
      end
 
 
 
            function y=f_ode(x1,x2,x3,choice,a,b)
 
            %  function 'f_ode' is the the functions 
            %  used in the LOTKA VOLTERA MODEL with RK, EULER, EULER-CAUCHY (above)
            %
            % <- Kalamaras Dimitris ->
 
 
            if nargin==0 error('Usage: f_ode(x1,x2,x3,choice,a,b ) , where xi real, choice {1 or 2} and a,b parameters'); end
 
            if choice==1
               y=a*x2 - x2*x3 ;  %% lotka VOLTERA MODEL f(x,y,z)=ay-yz %%
            elseif choice==2
               y=-b*x3+ x2*x3;  %% lotka voltera MODEL g(x,y,z)=-bz+yz %%
            else
               error('Usage: f_ode(x1,x2,x3,choice,a,b ) , where xi real, choice {1 or 2} and a,b parameters'); 
            end
 
 
 
 
      % Approximation of the integral of 4/(1+x^2) from zero to one via Runge - Kutta 4th class
      %
      % set y'=4/(1+x^2) and history y(0)=0, thus the integral is the solution of this ODE in 0<=x<=1 if y(1)=pi
      % where pi is the real value of the given integral. The method used in the integration of the ODE is Runge Kutta
      %
      % <- Kalamaras Dimitris ->
 
 
      clc;
      clear all;
      format long;
 
      h=input ('Enter step h for R-K : ');
      %h=0.01;
      tol=input('Enter desired precision TOL : ' );
 
      max=1/h;
      x=0; %
      y=0; % history zero;
      z=0;
      for i=1:max
 
         k1=h*f_rk(x,y);
      	k2=h*f_rk(x+h/2,y+k1/2);
      	k3=h*f_rk(x+h/2,y+k2/2);
      	k4=h*f_rk(x+h,y+k3);
         y=y+ (1/6) * (k1+2*k2+2*k3+k4);
         if abs(y-z)<= tol break; end
         z=y;
      	if x == 1
         disp('0k');break;
         end
      	x=i*h;
 
 
      end
      disp(' ');
      disp('RUNGE KUTTA RESULT:');
      sprintf(' y(1) =  %15.14f in %5d RUNGE KUTTA iterations', y, i-1)
 
      disp('while the real integral(4/(1+x^2)) is ' );
      syms r;
      real_value=int(4/(1+r^2),0,1)
 
 
      difference=abs(real_value-y)	
 
 
      x=0;
      y=0;
      z=0;
      for i=1:max
 
      	k1=h*f_rk(x,y);
      	k2=h*f_rk(x+h,y+k1);
         y=y+ (1/2)*(k1+k2);
         if abs(y-z)<= tol break; end
         z=y;
      	if x == 1
         disp('0k');break;
         end
         x=i*h;
 
 
      end
 
      disp(' ' );
      disp('EULER - HEUN RESULT:');
      sprintf(' y(1) = %15.14f in %5d EULER-HEUN iterations', y, i-1)
 
      difference=abs(y-real_value)	
 
 
 
            function z=f_rk(x,y)
 
            %  function 'f_rk' is the the function 
            %  used in the RungeKutta method above
            %
            z=4/(1+x^2);
 
 
 
 
 
 
      % Boundary Value D.E
      % Numerical Solution usinh the Shooting non-linear method with RUNGE KUTTA 4th class for
      % the solution of the initial value problems and Bisection and Newton Raphson method for t
      %
      %  y''=f(x,y,y') = 3y'-2y+x              % if you want another system, change f_ode1.m
      %  y(0)=1
      %  y(1)=2                   shooting: y'(0)=t						
      %
      %
      % <- Dimitris B. Kalamaras  ->
 
      clc;
      clear all;
 
      disp( 'SHOOTING METHOD with RUNGE-KUTTA 4th Class and NEWTON-RAPHSON or BISECTION ');
      disp( '  y''''=f(x,y,y'') = 3y''-2y+x');
      disp( '  y(0)=1,   y(1)=1    -->       shooting: y''(0)=t');
      disp( ' ');																			
 
      %DATA INPUT & INITIALIZATION
      h=input   ('RUNGE-KUTTA...: Enter step h for RK [0.001]: ');
      e=input   ('NEWTON-RAPHSON: Enter step e for NR [0.1]: ');
      t(1)=input('NR: TARGET....: Enter initial estimation of t : ');     % newton raphson only
      prec =input('NR & BISECTION: t approximation precision [<0.1]: ');
      a=input   ('BISECTION.... : Enter left interval margin a [0.1]: ');
      b=input   ('BISECTION ....: Enter right interval margin b [3]: ');
      max=input ('Enter max. amount of iterations for shooting [8] : ') ; 
      n=1/h;
      Target=2;
 
 
      %
      % Runge Kutta with Finite-Differences Newton-Raphson
      %
      for k=1:max
        % solve the equation P(t)=y(b;t)-B=0 with a method, f.ex. Newton Raphson to find a better t.
        % Thus we need the root of the P
        % I need the derivative Of P(t)=y(b,t)-B for Newton Raphson, thus i use dP/dt ~ ( P(t+e)-P(t) )/e  as an approxim.
        % for that reason we calculate first  y with t=t+e:
       x(1)=0;         %bored changing the damn zero-index matlab option...who cares whilst it works, ?
       y(1)=1;
       z(1)=t(k)+e;
 
      for i=1:n
          k1=h*f_ode1(x(i),y(i),z(i),1);
          l1=h*f_ode1(x(i),y(i),z(i),2);
          k2=h*f_ode1(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,1);
          l2=h*f_ode1(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,2);
          k3=h*f_ode1(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,1);
          l3=h*f_ode1(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,2);
          k4=h*f_ode1(x(i)+h, y(i)+k3, z(i)+l3,1);
          l4=h*f_ode1(x(i)+h, y(i)+k3, z(i)+l3,2);
 
          y(i+1)=y(i)+ (1/6)*(k1 + 2*k2 + 2*k3 + k4);
          z(i+1)=z(i)+ (1/6)*(l1 + 2*l2 + 2*l3 + l4);
 
          x(i+1)=x(i)+h;
 
      end
 
      Pte=(y(n+1)-Target);
 
       % Now calculate y with target t=t   
 
       x(1)=0;         %bored changing the damn zero-index matlab option...who cares whilst it works, ?
       y(1)=1;
       z(1)=t(k);
 
       for i=1:n
          k1=h*f_ode1(x(i),y(i),z(i),1);
          l1=h*f_ode1(x(i),y(i),z(i),2);
          k2=h*f_ode1(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,1);
          l2=h*f_ode1(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,2);
          k3=h*f_ode1(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,1);
          l3=h*f_ode1(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,2);
          k4=h*f_ode1(x(i)+h, y(i)+k3, z(i)+l3,1);
          l4=h*f_ode1(x(i)+h, y(i)+k3, z(i)+l3,2);
 
          y(i+1)=y(i)+ (1/6)*(k1 + 2*k2 + 2*k3 + k4);
          z(i+1)=z(i)+ (1/6)*(l1 + 2*l2 + 2*l3 + l4);
 
          x(i+1)=x(i)+h;
      end
       disp('SHOOTING with NEWTON-RAPHSON:');
       sprintf('In %2d repeat, for n= %6d  the desired y(n,t) has value %20.17f with t = %8.4f', k, n, y(n+1),t(k))
 
       Pt=y(n+1)-Target;  % function P(t) of the difference y(b,t)-target will be used in Newton Raphson to find a better t
        if abs(Pt) < prec 
         disp('Shooting Succeeded!');   t(k+1)=t(k); break;
        end;
 
      % approximation of the derivative
      dPt=(Pte-Pt)/e;
      t(k+1)=t(k)-Pt/dPt;
 
      figure;
      plot(x,y,'g:');
      title('SHOOTING & NR: Graph of y in between approximations of t ');
      end
      figure;
      plot(x,y,'r:');
      title('SHOOTING & NR: Final graph of y in the interval [0,1]');
      disp('NR Iterations Needed :' ); disp (k);
 
 
 
      disp('Wait...Writing plot for Newton-Raphson...');
      disp(' ');
 
 
 
 
 
 
 
      %
      % Runge Kutta with BISECTION
      %
 
      % FIRST CALCULATE y with target t=a , t=b, i.e. the left % right margins of the interval...
      % this is done because I need the value (sign) of the function Pt= y(b,t)-Target at t = a
      % and also the value at t=b
      for k=1:2
 
      	x(1)=0; 
      	y(1)=1;  
      	if k==1 	  z(1)=a;
      	  elseif k==2  z(1)=b;
      	end
 
      	for i=1:n
      	    k1=h*f_ode1(x(i),y(i),z(i),1);
      	    l1=h*f_ode1(x(i),y(i),z(i),2);
      	    k2=h*f_ode1(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,1);
      	    l2=h*f_ode1(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,2);
      	    k3=h*f_ode1(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,1);
      	    l3=h*f_ode1(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,2);
      	    k4=h*f_ode1(x(i)+h, y(i)+k3, z(i)+l3,1);
      	    l4=h*f_ode1(x(i)+h, y(i)+k3, z(i)+l3,2);
      	    y(i+1)=y(i)+ (1/6)*(k1 + 2*k2 + 2*k3 + k4);
      	    z(i+1)=z(i)+ (1/6)*(l1 + 2*l2 + 2*l3 + l4);
 
      	    x(i+1)=x(i)+h;
         end
 
      	if k==1           Pta=(y(n+1)-Target);   % Got it.
         	elseif k==2    Ptb=(y(n+1)-Target);   % Got it.
      	end;
 
      end
 
      t(1)=a+((b-a)/2^1)*sign(Pta)*sign(Pta); % calculate initial t with BISECTION
 
      % Now approximate t(k) with Bisection in max times, by calculating y, and testing the Pt=y(b,t)-Target function.
 
      for k=1:max
 
       % Now calculate y with target t=t(k) given from Bisection   
 
       x(1)=0;         %bored changing the damn zero-index matlab option...who cares whilst it works, ?
       y(1)=1;
       z(1)=t(k);
 
       for i=1:n
          k1=h*f_ode1(x(i),y(i),z(i),1);
          l1=h*f_ode1(x(i),y(i),z(i),2);
          k2=h*f_ode1(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,1);
          l2=h*f_ode1(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,2);
          k3=h*f_ode1(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,1);
          l3=h*f_ode1(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,2);
          k4=h*f_ode1(x(i)+h, y(i)+k3, z(i)+l3,1);
          l4=h*f_ode1(x(i)+h, y(i)+k3, z(i)+l3,2);
 
          y(i+1)=y(i)+ (1/6)*(k1 + 2*k2 + 2*k3 + k4);
          z(i+1)=z(i)+ (1/6)*(l1 + 2*l2 + 2*l3 + l4);
 
          x(i+1)=x(i)+h;
      end
       disp('Shooting with BISECTION :');
       sprintf('In %2d repeat, for n= %6d  the desired y(n,t) has value %20.17f with t = %8.4f', k, n, y(n+1),t(k))
 
       Ptk=y(n+1)-Target;  % function P(t(k)) of their difference will be used in BISECTION to find a better t
       if abs(Ptk)  < prec
         disp('Shooting Succeeded!');  t(k+1)=t(k);   break;
       end;
 
 
       % BISECTION ALGORITHM
 
       if Pta*Ptk >0 & Ptb*Ptk > 0 ...
             disp('No roots for P(t)=y(b,t)-Target in the given interval ');  break; 
       end; 
 
      t(k+1)=t(k)+((b-a)/2^(k+1))*sign(Pta)*sign(Ptk);
 
 
 
      figure;
      plot(x,y,'g:');
      title('SHOOTING & BISECTION: Graph of y in between approximations of t ');
      end
      figure;
      plot(x,y,'r:');
      title('SHOOTING & BISECTION: Final graph of y in the interval [0,1]');
      disp('BISECTION Iterations Needed :' ); disp (k);
 
 
 
 
       %for educational purposes plot y from 0 to 5
 
       x(1)=0;        y(1)=1;
       z(1)=t(k+1);
       h=0.001;
       n=5/h ;
       for i=1:n
          k1=h*f_ode1(x(i),y(i),z(i),1);
          l1=h*f_ode1(x(i),y(i),z(i),2);
          k2=h*f_ode1(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,1);
          l2=h*f_ode1(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,2);
          k3=h*f_ode1(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,1);
          l3=h*f_ode1(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,2);
          k4=h*f_ode1(x(i)+h, y(i)+k3, z(i)+l3,1);
          l4=h*f_ode1(x(i)+h, y(i)+k3, z(i)+l3,2);
 
          y(i+1)=y(i)+ (1/6)*(k1 + 2*k2 + 2*k3 + k4);
          z(i+1)=z(i)+ (1/6)*(l1 + 2*l2 + 2*l3 + l4);
 
          x(i+1)=x(i)+h;
 
      end
      figure;
      plot(x,y);
      title('BISECTION Final Plot of integral y of the given ODE for the interval [0,5] ');
 
 
 
            function y=f_ode1(x1,x2,x3,choice)
 
            %  function 'f_ode1' is the system of the functions F,G
            %  used in the with RK above
            %
            %  The initial valued ODE { y''=3y'-2y+x   y(0)=	1, y'(0)=t } 
            %  with y1=y, y2=y', y3=y'' becomes the system of 1st class ODEs: 
            %  -- { y1'=y2,  with y1(0)=1
            %  -- { y2'=y''=3y'-2y+x=3y2-2y+x   with y2(0)=t
            %
            % <- Kalamaras Dimitris ->
 
 
 
 
            if nargin==0 error('Usage: f_ode1(x1,x2,x3,choice ) , where xi real, choice {1 or 2} '); end
 
            if choice==1
               y=x3 ;  %%  %%
            elseif choice==2
               y=3*x3-2*x2+x1;  %% %%
            else
               error('Usage: f_ode1(x1,x2,x3,choice ) , where xi real, choice {1 or 2} '); 
            end
 
 
 
 
    * RANDOM NUMBER GENERATOR ADAPTED FROM NUMERICAL RECIPES
 
      function y=ran0(seed)
 
      % Returns a uniform random deviate between 0.0 and 1.0.  
      % The period of ran0 is 2^31~2.1x10^9
      %
      % Set or reset idum to any integer (except the unlikely value MASK) to initialize
      % the sequence; idum must not be altered between calls for successive deviates
      % in a sequence. Value 0 must not be allowed as seed and should not be expected as rand.num!
      % This is a peculiarity of all r.n.g of the form I(n+1)=a*I(n) (mod m)
      %
      % Taken from Numerical Recipes in C Internet Edition.
      % Adapted in Matlab for creating noise over a function.
      % View ran0.m for more details...if you have lost NR in C!
      %
      % <-  Kalamaras Dimitris ->
 
      format long;
      global idum;
 
      % Park & Miller "Minimal Standard" proposal
      if nargin==1 idum=seed; end
      idum;
      ia=16807;
      im=2147483647;  %2^31-1
      am=1.0/im;
      %MASK=123459876
            % if im=ia*iq+ir (Shcrage Method)   
             %  then ia*idum=a(z mod q)-r[z/q]  (+im if <0) 
              % for multiplying two 32-bit integers modulo a 32-bit constant,
               % without using any intermediates larger than 32 bits (including a sign bit)
 
      iq=127773;
      ir=2836;
                             % I should XOR with MASK to allow zero seed (read NR) but...
      %idum=xor(idum,MASK)    % unfortunately  XORing in matlab requires vectors x,y!
      k=round(idum/iq);        % Done with round (thus, less efficiently) to keep code NR-like
      idum=ia*(idum-iq*k)-ir*k; % It also seems that I am avoidind overflows...DO NOT USE MOD!
      if idum < 0                % If it is <0 im should be added.
         idum=idum+im;
      end
                                    %convert to fp
      y=am*idum;
      %idum=xor(idum,MASK)
 
 
    *
 
      % NUMERICAL ANALYSIS
      % This script integrates error function of a standard normal variable
      % manifestated in f_mc_erf.m  using the Mean Sample Monte Carlo method
      %
      % <- Kalamaras Dimitris ->
 
 
      clc;
      clear all;
      format long;
      disp('NUMERICAL INTEGRATION');
      disp('Monte Carlo Method of error function in (0,1)');
      tic;
      I_Estimator=zeros(5,1);
      %n=input('Enter Amount of Random Draws for Monte Carlo :');
      n=40000;
      trial=1;
      while  n < 40001
      sum=0;
      rand('state',0);
      for i=1:n
         x=rand;
         sum=sum+F_mc_erf(x);
      end
 
      I_Estimator(trial)=(1/n)*sum;
      sum=0;
      rand('state',0);
      for i=1:n
         x=rand;
         sum=sum+(F_mc_erf(x)-I_Estimator(trial))^2;
      end
 
      disp(sum);
 
      Var_Estimator=sqrt( (1/n) * (1/(n-1))*sum );
      Var_1=sqrt ( (1/n) * (1/(n-1)) * sum );	%??????????????????????? zero ????????
 
      disp(n);
      sprintf('The approximate value of the integral is I = %10.5f',I_Estimator(trial))
      sprintf('With Error of the I to be = %10.5f',Var_Estimator)
      sprintf('With Variance  of the I to be = %10.5f',Var_1)
      n=n+1000;
      trial=trial+1;
      end
 
 
      x=zeros(5,1);
 
      for i=1:5
         if i==1 x(i)=1000
         else x(i)=x(i-1)+1000;
         end
      end
 
      plot (x, I_Estimator, '.');
      toc;
 
 
 
 
 
 
 
      % MONTE CARLO INTEGRATION
      % This script integrates a given univariate function manifestated in f_mc1.m  
      % using the Monte Carlo method
      %
      % <- Kalamaras Dimitris ->
 
 
      clc;
      clear all;
      disp('NUMERICAL INTEGRATION');
      disp('Monte Carlo Method');
      n=input('Enter Amount of Random Draws for Monte Carlo :');
      x=zeros(n);
      sum=0;
 
      for i=1:n
         x(i)=rand;
         sum=sum+F_mc1(x(i));
      end
 
      I_Estimator=(1/n)*sum;
      sum=0;
      for i=1:n
         sum=sum+(F_mc1(x(i))-I_Estimator)^2;
      end
 
      Var_Estimator=(1/n)*(1/(n-1))*sum;
 
      sprintf('The approximate value of the integral is I = %10.5f',I_Estimator)
      sprintf('With Variance of the I to be = %10.5f',Var_Estimator)
 
 
 
            function y=f_mc1(x)
 
            %  function 'f_mc1' is the the host m file of the function which
            %  is being integrated by the "mc1" script
            %
            % <- Kalamaras Dimitris ->
 
 
            if nargin==0 x=0:.1:10; end
            %if nargin==0 error('Usage: f_mc1 ( x ) , where x whatever.'); end
            y=cos(x);
            %y=x+x^4;
 
 
 
 
 
 
 
 
      % This script integrates a given two-variate function manifestated in f_mc2.m  
      % using the Monte Carlo method
 
      % <- Kalamaras Dimitris ->
 
 
      clc;
      clear all;
      disp('NUMERICAL INTEGRATION');
      disp('Monte Carlo Method in 2 dimensions ');
      a1=input('Input left margin a1 for the first dimension :');
      b1=input('Input right margin b1 for the first dimension :');
      a2=input('Input left margin a2 for the second dimension :');
      b2=input('Input right margin b2 for the second dimension :');
      k=1/((b1-a1)*(b2-a2))
      n=input('Enter Amount of Random Draws for Monte Carlo :');
      %x=zeros(n);
      %y=zeros(n);
      sum=0;
 
      for i=1:n
      %   x(i)=(b1-a1)*rand-a1;        %transformers....enotheite!
      %   y(i)=(b2-a2)*rand-a2;
      x=(b1-a1)*rand-a1;
      y=(b2-a2)*rand-a2;
      %   sum=sum+F_mc2(x(i),y(i));
      	  sum=sum+F_mc2(x,y);
 
      end
 
      I_Estimator=(1/(k*n))*sum;
 
      %sum=0;
      %for i=1:n
      %   sum=sum+(F_mc1(x(i))-I_Estimator)^2;
      %end
 
      %Var_Estimator=(1/n)*(1/(n-1))*sum;
 
      sprintf('The approximate value of the integral is I = %10.5f',I_Estimator)
      %sprintf('With Variance of the I to be = %10.5f',Var_Estimator)
 
 
 
            function y=f_mc2(x,y)
 
            %  function 'f_mc2' is the the host m file of the 2dimensional function which
            %  is being integrated by the "mc2" script
            %
            % <- Kalamaras Dimitris ->
 
 
 
            if (nargin==0) x=0:.1:10; y=0:.1:10; end
            %if nargin==0 error('Usage: f_mc2 ( x,y ) , where x whatever.'); end
 
 
            %y=x^2+y^2*cos(x)-exp(x);
            y=x^2+y^2;