Hey there,

I've a numerical solved system of differential equations, which depend on one argument and one index. I can solve it, but when I try plot it I have this error: Error, (in plot) two lists or Vectors of numerical values expected.

Could anyone help me figure out what I'm doing wrong?

> restart;

> A := 115.1558549; B := .3050464658; n := 3; f0 := 0.5e-4;

>f:=theta->f0*(cos(arcsin(sin(theta)/n)))^2;

I0:=Ir(z)+sum(Is[k](z),k=1..20);

> alpha := [0, 1, 2, 3, 4, 5, 6];

Theta := [3*Pi*(1/180), 6*Pi*(1/180), 9*Pi*(1/180), 12*Pi*(1/180), 15*Pi*(1/180), 18*Pi*(1/180), 21*Pi*(1/180), 24*Pi*(1/180), 27*Pi*(1/180), 30*Pi*(1/180), 33*Pi*(1/180), 36*Pi*(1/180), 39*Pi*(1/180), 42*Pi*(1/180), 45*Pi*(1/180), 48*Pi*(1/180), 51*Pi*(1/180), 54*Pi*(1/180), 57*Pi*(1/180), 60*Pi*(1/180)];

>G:= theta->A*sin(theta)*cos(2*arcsin((sin(theta)/n)))/((1+sin(theta)^2/B^2)*cos(arcsin(sin(theta)/n)));

>for j from 1 to 7 do

d1 := diff(Ir(z), z) = -sum(G(Theta[k])*Ir(z)*Is[k](z)/I0,k=1..20)-alpha[j]*Ir(z)-sum(f(Theta[k])*Ir(z),k=1..20):

d2 := diff(Is[1](z), z) = G(Theta[1])*Ir(z)*Is[1](z)/I0-alpha[j]*Is[1](z)+f(Theta[1])*Ir(z):

d3 := diff(Is[2](z), z) = G(Theta[2])*Ir(z)*Is[2](z)/I0-alpha[j]*Is[2](z)+f(Theta[2])*Ir(z):

d4 := diff(Is[3](z), z) = G(Theta[3])*Ir(z)*Is[3](z)/I0-alpha[j]*Is[3](z)+f(Theta[3])*Ir(z):

d5 := diff(Is[4](z), z) = G(Theta[4])*Ir(z)*Is[4](z)/I0-alpha[j]*Is[4](z)+f(Theta[4])*Ir(z):

d6 := diff(Is[5](z), z) = G(Theta[5])*Ir(z)*Is[5](z)/I0-alpha[j]*Is[5](z)+f(Theta[5])*Ir(z):

d7 := diff(Is[6](z), z) = G(Theta[6])*Ir(z)*Is[6](z)/I0-alpha[j]*Is[6](z)+f(Theta[6])*Ir(z):

d8 := diff(Is[7](z), z) = G(Theta[7])*Ir(z)*Is[7](z)/I0-alpha[j]*Is[7](z)+f(Theta[7])*Ir(z):

d9 := diff(Is[8](z), z) = G(Theta[8])*Ir(z)*Is[8](z)/I0-alpha[j]*Is[8](z)+f(Theta[8])*Ir(z):

d10 := diff(Is[9](z), z) = G(Theta[9])*Ir(z)*Is[9](z)/I0-alpha[j]*Is[9](z)+f(Theta[9])*Ir(z):

d11 := diff(Is[10](z), z) = G(Theta[10])*Ir(z)*Is[10](z)/I0-alpha[j]*Is[10](z)+f(Theta[10])*Ir(z):

d12 := diff(Is[11](z), z) = G(Theta[11])*Ir(z)*Is[11](z)/I0-alpha[j]*Is[11](z)+f(Theta[11])*Ir(z):

d13 := diff(Is[12](z), z) = G(Theta[12])*Ir(z)*Is[12](z)/I0-alpha[j]*Is[12](z)+f(Theta[12])*Ir(z):

d14 := diff(Is[13](z), z) = G(Theta[13])*Ir(z)*Is[13](z)/I0-alpha[j]*Is[13](z)+f(Theta[13])*Ir(z):

d15 := diff(Is[14](z), z) = G(Theta[14])*Ir(z)*Is[14](z)/I0-alpha[j]*Is[14](z)+f(Theta[14])*Ir(z):

d16 := diff(Is[15](z), z) = G(Theta[15])*Ir(z)*Is[15](z)/I0-alpha[j]*Is[15](z)+f(Theta[15])*Ir(z):

d17 := diff(Is[16](z), z) = G(Theta[16])*Ir(z)*Is[16](z)/I0-alpha[j]*Is[16](z)+f(Theta[16])*Ir(z):

d18 := diff(Is[17](z), z) = G(Theta[17])*Ir(z)*Is[17](z)/I0-alpha[j]*Is[17](z)+f(Theta[17])*Ir(z):

d19 := diff(Is[18](z), z) = G(Theta[18])*Ir(z)*Is[18](z)/I0-alpha[j]*Is[18](z)+f(Theta[18])*Ir(z):

d20 := diff(Is[19](z), z) = G(Theta[19])*Ir(z)*Is[19](z)/I0-alpha[j]*Is[19](z)+f(Theta[19])*Ir(z):

d21 := diff(Is[20](z), z) = G(Theta[20])*Ir(z)*Is[20](z)/I0-alpha[j]*Is[20](z)+f(Theta[20])*Ir(z):

dsys := {d1, d10, d11, d12, d13, d14, d15, d16, d17, d18, d19, d2, d20, d21, d3, d4, d5, d6, d7, d8, d9}:

dSol[j] := dsolve({op(dsys), Ir(0) = 1, Is[1](0) = 0.1e-1, Is[2](0) = 0.1e-1, Is[3](0) = 0.1e-1, Is[4](0) = 0.1e-1, Is[5](0) = 0.1e-1, Is[6](0) = 0.1e-1, Is[7](0) = 0.1e-1, Is[8](0) = 0.1e-1, Is[9](0) = 0.1e-1, Is[10](0) = 0.1e-1, Is[11](0) = 0.1e-1, Is[12](0) = 0.1e-1, Is[13](0) = 0.1e-1, Is[14](0) = 0.1e-1, Is[15](0) = 0.1e-1, Is[16](0) = 0.1e-1, Is[17](0) = 0.1e-1, Is[18](0) = 0.1e-1, Is[19](0) = 0.1e-1, Is[20](0) = 0.1e-1}, [Ir(z), Is[1](z), Is[2](z), Is[3](z), Is[4](z), Is[5](z), Is[6](z), Is[7](z), Is[8](z), Is[9](z), Is[10](z), Is[11](z), Is[12](z), Is[13](z), Is[14](z), Is[15](z), Is[16](z), Is[17](z), Is[18](z), Is[19](z), Is[20](z)], numeric);

end do:

>for j from 1 to 7 do

dSol[j](0.4);

as:='as':

for l from 1 to 20 do

as[l]:=[Theta[l],rhs(dSol[j](0.4)[2+l])];

od:

plo[j]:=convert(as,listlist);

od:

>plot(plo[2],plo[1]);

Error, (in plot) two lists or Vectors of numerical values expected