graph did not plot, how to plot ?

function vdpd_qo_syn

clear all

clear all

figure(1)

%Computational Length and step size

n=100; %total data legnth (n/h)

h=0.01;

%%%%%%%%%%%%%%%%%%%system 1%%%%%%%%%%%%%%%%

alpha= 0.01;

beta = 9.74;

omegaf= 0.3;

gamma= 0.2;

W1= 0.3;

ff1 =0.4999;

%% Initial condition

%x10 = 0.0; x20 = 0.0;

x10 = 0.11; x20 = 0.2;

%%%%%%%%%%%%%%%system 2%%%%%%%%%%%%%%%%

e = 0.22;

a= 0.956;

b = 0.149;

ff2 =.59;

W2= 0.3;

%% Initial condition

%x30 = 0.0; x40 = 0.0;

x30 = 0.11; x40 = 0.2;

%%%%%%%%%%%%%%%%%%coupling %%%%%%%%%%%%%%%%

%cou=0.2;

%e1 = x3-x1; e2 = x4-x2;

%time step and initial condition

tspan = 0:h:n;

y0 = [x10; x20; x30; x40];

%[t,y] = ode45(@(t,x) f(t,x,a,b,beta,rho,ff1,W1,W2,ff2,cou),tspan,y0);

%[t,y] = ode45(@(t,x) f(t,x,alpha,omegaf,beta,gamma,ff1,ff2,W1,W2,epsilon,alpha1,beta1,cou),tspan,y0);

[t,y] = ode45(@(t,x) f(t,x,alpha,beta,omegaf,gamma,ff1,W1,e,a,b,ff2,W2),tspan,y0);

%[t,y] = ode45(@(t,x) f(t,x,alpha,beta,omegaf,gamma,ff1,W1,e,a,b,ff2,W2,cou),tspan,y0);

%% transient%%%%%%%%%%%%%%%%remove%%%%%%%%%%%%%%%

k=n/h;

k1=(k/2:k);

size(k1);

x1=y(k1,1); x2=y(k1,2); x3=y(k1,3); x4=y(k1,4); ti=t(k1);

%plot the variable

%plot(x1,x2)

subplot(2,2,1); plot(x1,x2); xlabel('x1'); ylabel('x2'); title('Quintic oscillator ');

subplot(2,2,2); plot(x3,x4); xlabel('y1'); ylabel('y2'); title('VdpD oscillator');

subplot(2,2,3); plot(x1,x3); xlabel('x1'); ylabel('y1'); title('Synchronization');

subplot(2,2,4); plot(x2,x4); xlabel('x2'); ylabel('y2'); title('Synchronization');

figure(2)

subplot(2,2,3);plot(ti,x1,'r',ti,x3,'b'); xlabel('Time(sec)'); ylabel('x1,y1');

subplot(2,2,4);plot(ti,x2,'r',ti,x4,'b'); xlabel('Time(sec)'); ylabel('x2,y2');

subplot(2,2,1);plot(ti,x1,'r'); xlabel('Time(sec)'); ylabel('x1');

subplot(2,2,2);plot(ti,x2,'b'); xlabel('Time(sec)'); ylabel('y1');

fprintf('Total length %d and the code taken transient after %d', k,k/2)

function dy = f(t,y,alpha,beta,omegaf,gamma,ff1,W1,e,a,b,ff2,W2)

%function dy = f(t,y,alpha,beta,omegaf,gamma,ff1,W1,e,a,b,ff2,W2,cou)

x1 = y(1); x2 = y(2); x3 = y(3); x4 = y(4);

%sys-I

dx1=x2;

dx2=-alpha*x2-omegaf*x1-beta*x1.^3- gamma*x1.^5+ff1*cos(W1*t);

%Sys II

u1 = -x3+x1-x4+x2;

u2 = -(omegaf-alpha +beta*x1.^2 +gamma*x1.^4)*x1 -b*x2.^3-(e+alpha)*x3 +e*x2.^2*x4 -ff2*sin(W2*t)-ff1*cos(W1*t);

dx3=x4+u1;

dx4=e*(1-x3.^2)*x4-a*x3+b*x3.^2+ff2*cos(W2*t)+u2;

dy = [dx1; dx2; dx3; dx4];

Answer 1

[Voss]

When I call the function, it plots:

% call the function vdpd_qo_syn:
vdpd_qo_syn()

Total length 10000 and the code taken transient after 5000

% define the function vdpd_qo_syn:
function vdpd_qo_syn
clear all
clear all
figure(1)
%Computational Length and step size
n=100;   %total data legnth (n/h)
h=0.01;
%%%%%%%%%%%%%%%%%%%system 1%%%%%%%%%%%%%%%%   
       
        alpha= 0.01;
         beta = 9.74;
        omegaf= 0.3;
        gamma= 0.2;
        W1= 0.3;
        ff1 =0.4999;
      %% Initial condition
      %x10 = 0.0; x20 = 0.0;
      x10 = 0.11; x20 = 0.2; 
    %%%%%%%%%%%%%%%system 2%%%%%%%%%%%%%%%%   
       e = 0.22;
       a= 0.956;
       b = 0.149;
       ff2 =.59;
       W2= 0.3;
         
      %% Initial condition
       %x30 = 0.0; x40 = 0.0; 
     x30 = 0.11; x40 = 0.2; 
     
  %%%%%%%%%%%%%%%%%%coupling %%%%%%%%%%%%%%%%   
  %cou=0.2;
  %e1 = x3-x1; e2 = x4-x2;
  
  %time step and initial condition
      tspan = 0:h:n;  
       y0 = [x10; x20; x30; x40];
       %[t,y] = ode45(@(t,x) f(t,x,a,b,beta,rho,ff1,W1,W2,ff2,cou),tspan,y0);   
      %[t,y] = ode45(@(t,x) f(t,x,alpha,omegaf,beta,gamma,ff1,ff2,W1,W2,epsilon,alpha1,beta1,cou),tspan,y0);   
 [t,y] = ode45(@(t,x) f(t,x,alpha,beta,omegaf,gamma,ff1,W1,e,a,b,ff2,W2),tspan,y0);   
 %[t,y] = ode45(@(t,x) f(t,x,alpha,beta,omegaf,gamma,ff1,W1,e,a,b,ff2,W2,cou),tspan,y0);   
 
      
%% transient%%%%%%%%%%%%%%%%remove%%%%%%%%%%%%%%%
  k=n/h;
  k1=(k/2:k);
  size(k1);
  x1=y(k1,1); x2=y(k1,2);  x3=y(k1,3); x4=y(k1,4); ti=t(k1);
  
   %plot the variable
   %plot(x1,x2)
   
   subplot(2,2,1); plot(x1,x2); xlabel('x1'); ylabel('x2'); title('Quintic oscillator '); 
   subplot(2,2,2); plot(x3,x4); xlabel('y1'); ylabel('y2'); title('VdpD oscillator'); 
   subplot(2,2,3); plot(x1,x3); xlabel('x1'); ylabel('y1'); title('Synchronization');
   subplot(2,2,4); plot(x2,x4); xlabel('x2'); ylabel('y2'); title('Synchronization'); 
   
   figure(2)
   subplot(2,2,3);plot(ti,x1,'r',ti,x3,'b'); xlabel('Time(sec)'); ylabel('x1,y1'); 
   subplot(2,2,4);plot(ti,x2,'r',ti,x4,'b'); xlabel('Time(sec)'); ylabel('x2,y2');
   subplot(2,2,1);plot(ti,x1,'r'); xlabel('Time(sec)'); ylabel('x1'); 
   subplot(2,2,2);plot(ti,x2,'b'); xlabel('Time(sec)'); ylabel('y1');
   
   fprintf('Total length %d and the code taken transient after %d', k,k/2)
end


%define the function f:
function dy = f(t,y,alpha,beta,omegaf,gamma,ff1,W1,e,a,b,ff2,W2) 
   %function dy = f(t,y,alpha,beta,omegaf,gamma,ff1,W1,e,a,b,ff2,W2,cou)        
   x1 = y(1);    x2 = y(2);   x3 = y(3);    x4 = y(4);
   
 
  
   %sys-I
   	
    dx1=x2;
	dx2=-alpha*x2-omegaf*x1-beta*x1.^3- gamma*x1.^5+ff1*cos(W1*t);
  %Sys II
    u1 = -x3+x1-x4+x2;
    u2 = -(omegaf-alpha +beta*x1.^2 +gamma*x1.^4)*x1 -b*x2.^3-(e+alpha)*x3 +e*x2.^2*x4 -ff2*sin(W2*t)-ff1*cos(W1*t);
          dx3=x4+u1;
    dx4=e*(1-x3.^2)*x4-a*x3+b*x3.^2+ff2*cos(W2*t)+u2;
     dy = [dx1; dx2; dx3; dx4];
end