Plot multiple solutions of a system of differential equations with different values of constant c by using NDsolve

I am traying to solve system of differential equations with different values of constant c then plot the solution with different values of the constant c. where c= 2.2758, 2.04822, 1.82064, 1.59306, 1.36548 So what i want is plot a variety of different graphs over the range of c and if it possible to get each graphs with different color.

my code

Clear[t]
τ = 13.8;
ω0 = 1;
r = 0.7071;
ℏ = 1.05457173*10^-34;
k = 1666666.667;
z = 0;
Ω = 2.2758;
τ = 13.8;
T2 = 200;
Δ = 1.7758;
ω = 0.5;
Table[c {c, {2.2758, 2.04822, 1.82064, 1.59306, 1.36548}}];

system1 = {u'[t] == Δ*v[t],
v'[t] == -Δ*u[t] -c*E^(-(r^2/ω0^2) - ((t^2*(1.177^2))/τ^2))*w[t],
w'[t] == c*E^(-(r^2/ω0^2) - ((t^2*(1.177^2))/τ^2))*v[t]};
initialvalues1 = {u[-20] ==  0, v[-20] ==  0, w[-20] == -  1};
sol1 = NDSolve[Join[system1, initialvalues1], {u[t], v[t], w[t]}, {t, -20, 20}];

F = (-10^33*c*ℏ*r*u[t])/ω0^2*E^(-(r^2/ω0^2) - ((t^2*(1.177^2))/τ^2));
P = Plot[Evaluate[F /. sol1], {t, -20, 20}, FrameLabel -> {"t(fs)", "FT"},
Frame -> True, PlotRange -> All, FrameTicks -> All]

• Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. – bbgodfrey Dec 26 '14 at 23:06
• All your questions seem strongly related, but you've been posting under other user handle before: mathematica.stackexchange.com/users/9850/saysics – Dr. belisarius Dec 27 '14 at 0:33
• Under that old handle you never voted and never accepted an answer. I hope you'll become a better citizen now. – Dr. belisarius Dec 27 '14 at 0:37

It appears that much of what you've posted doesn't work in its present form but it can be simplified greatly in a working form. I've highlighted the key changes that need to be made (with some explanation) and then have put it all together at the end.

Highlights

I'm not sure what you were trying to do using Table to generate the c values, but it seems that the simplest way would be to just include the values you want in a list like so:

cVals = {2.2758, 2.04822, 1.82064, 1.59306, 1.36548};


The next step is to solve the system of differential equations for different values of the parameter c. Since each of these functions u,v, and w are dependent on t, this is a good task for Mathematica's built in function ParametricNDSolveValue. There's a lot of great information in the official documentation for doing the type of thing you want to do. The key here is that you define c to be the parameter that varies between solutions.

pfun = ParametricNDSolveValue[Join[system1,initialvalues1],{u,v,w},{t,-20,20},{c}];


At the end of the day, it appears that the function F is really what is needed. We can define it as a function of c as follows.

F[t_,c_]:=With[{val = First@pfun[c][t]}, (-10^33*c*ℏ*r*val)/ω0^2*E^(-(r^2/ω0^2) - ((t^2*(1.177^2))/τ^2))];�������������������


The val=First@pfun[c][t] tells Mathematica that you want to evaluate the function u for a certain value of c and t. If you wanted to instead evaluate the value of v or w you would change First@ to Second@ or Third@ respectively because that's the order of variables given in pfun above.

Finally, the solutions for the desired values of c can be plotted. Here is where you want to use Table to plot solutions for each of the c values.

Plot[Evaluate[Table[F[t, c], {c, cVals}]], {t, -20, 20},PlotRange -> All, Frame -> True, FrameTicks -> All,FrameLabel -> {"t(fs)", "FT"}]�


Final Solution

All that's left now is to put it all together!

Clear[t]
(* define some constants *)
τ = 13.8;
ω0 = 1;
r = 0.7071;
ℏ = 1.05457173*10^-34;
k = 1666666.667;
z = 0;
Ω = 2.2758;
τ = 13.8;
T2 = 200;
Δ = 1.7758;
(* define the values of c *)
cVals = {2.2758, 2.04822, 1.82064, 1.59306, 1.36548};
(* define system of equations *)
system1 = {u'[t] == Δ*v[t], v'[t] == -Δ*u[t] -c*E^(-(r^2/ω0^2) - ((t^2*(1.177^2))/τ^2))*w[t], w'[t] ==
c*E^(-(r^2/ω0^2) - ((t^2*(1.177^2))/τ^2))*v[t]};
initialvalues1 = {u[-20] == 0, v[-20] == 0, w[-20] == -1};
(* compute the solution *)
pfun = ParametricNDSolveValue[Join[system1, initialvalues1], {u, v, w}, {t,-20, 20}, {c}];
F[t_, c_] := With[{val = First@pfun[c][t]}, (-10^33*c*ℏ*r*val)/ω0^2*E^(-(r^2/ω0^2) - ((t^2*(1.177^2))/τ^2))];
(* plot the solutions for the different c values *)
Plot[Evaluate[Table[F[t, c], {c, cVals}]], {t, -20, 20}, PlotRange -> All, Frame -> True, FrameTicks -> All, FrameLabel -> {"t(fs)", "FT"}]
(* celebrate! *)


You can also get a legend for the plot with the optional Plot argument:
PlotLegends->Map[ToString,cVals]