PLOT IS NOT SHOWING ANYTHING

Hello everyone,

I have a question regarding nothing showing up in my graph. How can I fix it?

My code is in the following:

clear all
clc
format long
%----parameters----%
global T0 Rg E1 E2 E_2 Eg k10 k20 k2_0 kgo Deo k1 k2 k_2 Ke Kg Cb Ct Ce a
a=1;
T0=800;%k Temperature
Rg=8.314;%J/k*mol; 
E1=10^5;%J/mol;
E2=10^5;%J/mol;
E_2=10^2;%J/mol;
Eg=10^4;%J/mol
k10=0.1;
k20=0.1;
k2_0=0.1;
kgo=0.1;
Deo=10^-10;
k1 = k10*exp(-E1/(Rg*T0));
k2 = k20*exp(-E2/(Rg*T0));
k_2 = k2_0*exp(-E_2/(Rg*T0));
Ke=(2.01*10^(-6))*T0^2-(1.43*10^(-3))*T0+0.2544;
Kg=kgo*exp(Eg/(Rg*T0));
Cb =0; 
Ct=20;
Ce=(7.67131719*10^(-47))*T0^(14.5144484);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r=linspace(0,100e-4,5);
t=linspace(0,300,100);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n = numel(r);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CC0 = zeros(n,1); 
CO0 = ones(n,1);
CD0 = ones(n,1); 
CM0 = zeros(n,1);
CA0 = zeros(n,1); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CO0(1)=15;
CD0(1)=5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y0 = [CC0;CO0;CD0;CM0;CA0];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[T, Y] = ode15s(@(t,y)fun(t,y,r,n),t,y0)%T is time
%%%%%%%%PLOT%%%%%%%%%%%%%%%%%%%
Y=length(r);
  plot(r,Y)
  set(gca,'YLim',[0 1e-2],'xLim',[0 5e-2]);
  xlabel('Radius')
  ylabel('concentration')
  legend('CC','CO','CD','CM','CA')

Function:

function DyDt=fun(t,y,r,n)
global T0 Rg E1 E2 E_2 Eg k10 k20 k2_0 kgo Deo k1 k2 k_2 Ke Kg Cb Ct Ce a
 
CC = zeros(n,1);
CO = zeros(n,1);
CD = zeros(n,1);
CM = zeros(n,1);
CA = zeros(n,1);
DCCDt = zeros(n,1);
DCODt = zeros(n,1);
DCDDt = zeros(n,1);
DCMDt = zeros(n,1);
DCADt = zeros(n,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
DyDt = zeros(5*n,1);   
rhalf = zeros(n-1,1); 
DCCDr = zeros(n-1,1);  
D2CCDr2 = zeros(n-1,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CC = y(1:n);
CO = y(n+1:2*n);
CD = y(2*n+1:3*n);
CM = y(3*n+1:4*n);
CA = y(4*n+1:5*n);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rhalf(1:n-1)=(r(1:n-1)+r(2:n))/2; 
%%%%%%%B.C.Left%%%%%%%%%%%%%%%%%
DCCDr(1)=0;
D2CCDr2(1) = 1.0/(rhalf(1)-r(1))*(CC(2)-CC(1))/(r(2)-r(1));
%%%%%%%B.C.Right%%%%%%%%%%%%%%%%%
DCCDt(n)=Cb
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
DCODt(n) =-k1*(CO(n)^2 -( 1/Ke)*CD(n)^2);
DCDDt(n) =k1*(CO(n)^2 - (1/Ke)*CD(n)^2)-k2*CD(n)+k_2*CA(n)*CM(n)*CC(n);
DCMDt(n) =k2*CD(n)-k_2*CA(n)*CM(n)*CC(n);
DCADt(n) =k2*CD(n)-k_2*CA(n)*CM(n)*CC(n);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=2:n-1
  DCCDr(i) = ((r(i)-r(i-1))/(r(i+1)-r(i))*(CC(i+1)-CC(i))+(r(i+1)-r(i))/(r(i)-r(i-1))*(r(i)-CC(i-1)))/(r(i+1)-r(i-1));
  D2CCDr2(i) = (rhalf(i)*(CC(i+1)-CC(i))/(r(i+1)-r(i))-rhalf(i-1)*(CC(i)-CC(i-1))/(r(i)-r(i-1)))/(rhalf(i)-rhalf(i-1));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:n-1
    
  De=Deo*exp(a*CM(i)/Ct);
  DCCDt(i) =2*De*DCCDr(i)+De*r(i)*D2CCDr2(i)+k2*CD(i)-k_2*CA(i)*CM(i)*CC(i);
  DCODt(i) =-k1*(CO(i)^2 -( 1/Ke)*CD(i)^2);
  DCDDt(i) =k1*(CO(i)^2 - (1/Ke)*CD(i)^2)-k2*CD(i)+k_2*CA(i)*CM(i)*CC(i);
  DCMDt(i) =k2*CD(i)-k_2*CA(i)*CM(i)*CC(i);
  DCADt(i) =k2*CD(i)-k_2*CA(i)*CM(i)*CC(i);
  
end
DyDt = [DCCDt;DCODt;DCDDt;DCMDt;DCADt]; 
end
 

Can anyone help? Thank you very much!

[darova]

Do you have any idea? Where is a mistake can be?

[S H]

I think the mistake is

remove Y=length(r)

and plot(t,Y) instead of plot(r, Y)

but I would like to plot the concentration of different species vs. time in separate graphs and the concentration of different species vs radius in in separate graphs.

I have no idea how to do it? Do you have any idea?

Thanks!

[Rik]

What exactly do you want to use as x and what as y?

[S H]

Hello Rik,

The picture above is what I would like to plot.

Answer 1

[darova]

Use subplot to draw curvees in separate windows

subplot(2,1,1)
plot(T,Y(:,1),'r')
subplot(2,1,2)]
plot(T,Y(:,2),'r')

[S H]

Hello darova,

I used the method you suggested. "plot(T, Y(:,1),'r')" is the conc. vs. time for one species at a fixed radius. I was wondering if we can plot the conc. vs. time for one species at different radius in the same graph. Also, the conc. vs. radius for one species at different times in the same graph. The picture below is what exactly I would like to plot. Thank you very much for your help!

[Rik]

Did you do a simple Google search for how to plot multiple graphs in a single axes? That should turn up the hold function.

[darova]

BUt your radius is not changing. I don't understand this graph: time and conc are changing (r variable is fixed in your code)