plot from ode45

function kk1

clear all

close all

options = odeset('RelTol',1e-9,'AbsTol', 1e-16);

a0=0.1;

b0=0.1;

p=0.2;

theta=30:70;

alpha=sqrt(1-0.25*p^2);

a=(alpha-sqrt(1-(p)^2))/sqrt(alpha^2-0.25*sin(theta));

sol =ode45(@q,[0 1200],[0.01 0.01 0.01 0.01 0.01 0], options)

end

% output

x1=sol(1).y(1,:);

y1=sol(1).y(2,:);

plot(x,y)

function dy=q(x)

f=sqrt(1+2*x);

b=(1+2*x)^-1.5;

m=16*alpha^2-0.50*sin(theta);

m=4*(2*alpha^2-8*(alpha*sqrt(1-p^2)));

o=16*sqrt(alpha^2-0.50*sin(theta));

A=y*(1+2.*x)^-1.5.*(2+a);

B=-1i*4*y*(0.633)^2*(alpha.^2-0.25*sin(theta))+2*sqrt(alpha^.2-0.25.*sin(theta))*(alpha-sqrt(1-(p).^2));

C=y*0.633^2*(m+n+o).*(alpha-2*sqrt(1-(p).^2));

D=- y*((2+a)^-1)*f^-2*(1.5+1i*(9/8)*f*b);

E=-0.25*0.633*0.633*y*(4-5*p^2+4*(1-exp(0.18/4)));

F=(1/8)*p^2*b0*a0*((1.70*sin(theta)+3.5*sqrt(alpha^2-0.25*sin(theta))+2*(alpha-sqrt(1-(p)^2))));

dy=A+B+C+D+E+F;

end

end

pl you are req to plot x vs y by varing theta from 30 to 70

Answer 1

[Yongjian Feng]

Too many errors in your code. I tried to clean some, but it still has trouble to figure out what y is:

function  kk1
    clear all 
    close all
    options = odeset('RelTol',1e-9,'AbsTol', 1e-16);

    a0=0.1;
    b0=0.1;
    p=0.2;
    theta=30:70;
    alpha=sqrt(1-0.25*p^2);
    a=(alpha-sqrt(1-(p)^2))./sqrt(alpha^2-0.25*sin(theta));

    sol =ode45(@q,[0 1200],[0.01 0.01 0.01 0.01 0.01 0], options)
    x1=sol(1).y(1,:); 
    y1=sol(1).y(2,:);
    plot(x,y)
end
% % output
% x1=sol(1).y(1,:); 
% y1=sol(1).y(2,:);
% plot(x,y) 
%  
function dy=q(x, alpha)
    f=sqrt(1+2*x);
    b=(1+2*x)^-1.5;
    theta=30:70;
    p=0.2;
    m=16*alpha.^2-0.50*sin(theta);
    m=4*(2*alpha.^2-8*(alpha*sqrt(1-p^2)));
    o=16*sqrt(alpha.^2-0.50*sin(theta));
    A=y*(1+2.*x)^-1.5.*(2+a); %What is y here?
    B=-1i*4*y*(0.633)^2*(alpha.^2-0.25*sin(theta))+2*sqrt(alpha^.2-0.25.*sin(theta))*(alpha-sqrt(1-(p).^2));
    C=y*0.633^2*(m+n+o).*(alpha-2*sqrt(1-(p).^2));
    D=- y*((2+a)^-1)*f^-2*(1.5+1i*(9/8)*f*b);
    E=-0.25*0.633*0.633*y*(4-5*p^2+4*(1-exp(0.18/4)));
    F=(1/8)*p^2*b0*a0*((1.70*sin(theta)+3.5*sqrt(alpha^2-0.25*sin(theta))+2*(alpha-sqrt(1-(p)^2))));
    dy=A+B+C+D+E+F;
end
% end

[shiv gaur]

dy is the out put from that can be find out from ode 45

here dy is differential

[shiv gaur]

edit code

function kk1

clear all

close all

options = odeset('RelTol',1e-9,'AbsTol', 1e-16);

a0=0.1;

b0=0.1;

p=0.2;

theta=30:70;

alpha=sqrt(1-0.25*p^2);

a=(alpha-sqrt(1-(p)^2))./sqrt(alpha^2-0.25*sin(theta));

sol =ode45(@q,[0 1200],0.01] options)

x1=sol(1).y(1,:);

y1=sol(1).y(2,:);

plot(x,y)

end

% % output

% x1=sol(1).y(1,:);

% y1=sol(1).y(2,:);

% plot(x,y)

%

function dy=q(x, alpha)

f=sqrt(1+2*x);

b=(1+2*x)^-1.5;

theta=30:70;

p=0.2;

m=16*alpha.^2-0.50*sin(theta);

m=4*(2*alpha.^2-8*(alpha*sqrt(1-p^2)));

o=16*sqrt(alpha.^2-0.50*sin(theta));

A=y*(1+2.*x)^-1.5.*(2+a); %What is y here?

B=-1i*4*y*(0.633)^2*(alpha.^2-0.25*sin(theta))+2*sqrt(alpha^.2-0.25.*sin(theta))*(alpha-sqrt(1-(p).^2));

C=y*0.633^2*(m+n+o).*(alpha-2*sqrt(1-(p).^2));

D=- y*((2+a)^-1)*f^-2*(1.5+1i*(9/8)*f*b);

E=-0.25*0.633*0.633*y*(4-5*p^2+4*(1-exp(0.18/4)));

F=(1/8)*p^2*b0*a0*((1.70*sin(theta)+3.5*sqrt(alpha^2-0.25*sin(theta))+2*(alpha-sqrt(1-(p)^2))));

dy=A+B+C+D+E+F;

end

% end