%        Compare  euler_multidim.m and rk4_multidim 
%            with the MATLAB provided packages
%                    ode45 and ode23
%        in the case of the anharmonic hoscillator
%
%                  d^2 x/dt^2=-x-x^3

clc                                   % clean screen
clear all                             % clear previously defined variables


% Parameters
 tinit=0.0;                            % initial time
 xinit=[0.2; 2.3];                     % initial position as raw vector
 tfin=2*pi;                            % final time
 steps=1000;                           % number of steps
 dim=2;                                % spatial dimensions

% define a rectangular matrix to store results
 y=zeros(dim+1,steps+1,2);    

% define input for euler_multidim (as raw vector)
 fun = @(t,y1,y2)  [0 1 ; -1 0]*[y1+y1^3;y2];
  
% define input for the MATLAB packages (as raw vector)
% fun_matlab = @(t,y)  [y(2) ; -y(1)-y(1)^3];                                          

% Solve with euler_multidim
 [y(1,:,1),y(2:dim+1,:,1)]=euler_multidim(fun,dim,tinit,xinit,tfin,steps);
 euler_steps=length(y(1,:,1));
 
% Solve with rk4_multidim
 [y(1,:,2),y(2:dim+1,:,2)]=rk4_multidim(fun,dim,tinit,xinit,tfin,steps);
 rk4_steps=length(y(1,:,2));
   
% Solve with ode45
  y45=[];                            % initializes a new vector
  [t,x]=ode45(@anharmonic, [tinit tfin],xinit);
  y45(1,:)=t;
  y45(2:dim+1,:)=x';                 % Output with time evolution in rows. 
  ode45_steps=length(y45(1,:));
  
% Solve with ode23
  y23=[];                            % initializes a new vector
  [t,x]=ode45(@anharmonic, [tinit tfin],xinit);
  y23(1,:)=t;
  y23(2:dim+1,:)=x';                 % Output with time evolution in rows.
  ode23_steps=length(y23(1,:));                                                     

%%% Plotting

% opens a new figure window 
 figure(1)         
 subplot(121)                        % divides window into a 1x2 matrix 
                                     % and assigns next plot to the 
                                     % left-hand frame (frame 1)
                                    
 plot(y(dim,:,1),y(dim+1,:,1))       % plot second versus first coordinates
 xlabel('x1')                        % label the x-axis
 ylabel('x2')                        % label the y-axis
 mesh_size1=['Mesh size:' sprintf('%0.3f',1/euler_steps)];
 text(floor(mean(y(dim,:,1))) ,floor(mean(y(dim+1,:,1))), mesh_size1 )
 title('Euler phase plane plot')     % gives a title to the plot
   
 subplot(122)                        % assigns next plot to right-hand 
                                     % frame (frame 2)
   
   % Plot the first coordinate versus time
   
 plot(y45(dim,:),y45(dim+1,:));
 xlabel('x1')                        % label the x-axis
 ylabel('x2')                        % label the y-axis
 mesh_size2=['Mesh size:' sprintf('%0.3f',1/ode45_steps)];
 text(-1 ,0, mesh_size2 )
 title('Ode45 phase plane plot')     % gives a title to the plot
   
 subplot(111)                        % return plotting mode to one 
 
% opens a new figure window
 figure(2)                
 subplot(121)                        % divides window into a 1x2 matrix 
                                     % and assigns 1st plot to the 
                                     % left-hand frame (frame 1)
                                    
 plot(y(dim,:,2),y(dim+1,:,2))       % plots second versus first coordinate
 xlabel('x1')                        % label the x-axis
 ylabel('x2')                        % label the y-axis
 mesh_size3=['Mesh size:' sprintf('%0.3f',1/rk4_steps)];
 text(floor(mean(y(dim,:,2))) ,floor(mean(y(dim+1,:,2))), mesh_size3 )
 title('Runge-Kutta phase plane plot')     % gives a title to the plot
   
 subplot(122)                        % assigns 2nd plot to the right-hand 
                                     % frame (frame 2)
                                     
   
% Plot the first coordinate versus time
   
 plot(y23(dim,:),y23(dim+1,:))  % plots second versus first coordinate
 xlabel('x1')                        % label the x-axis
 ylabel('x2')                        % label the y-axis 
 mesh_size4=['Mesh size:' sprintf('%0.3f',1/ode23_steps)];
 text(-1 ,0, mesh_size4 )
 title('Ode23 phase plane plot')     % gives a title to the plot

 