I need to solve two state variables using ODEs. I have set up the problem using NDSolve and proper boundary conditions. I am getting errors while Mathematica is trying to solve the set of 2 ODEs. I would appreciate any help in figuring out the problem.
R = 518.3; k = 1.304; tp = 10.; z =.;
V1 = .34; p1 = 25 10^6; T1 = 300;
\[Mu] = .002; Lpipe = 10.; dpipe = 3*25.4 10^-3; Rf = (
128 \[Mu] Lpipe)/(\[Pi] dpipe^4);
V = V1; \[Rho]a = 1.225; Kf = 1.8 10^9;
c = Sqrt[Kf/\[Rho]f]; \[Rho]f = 872; Ca = V/(\[Rho]a c^2);
Lt = 3*25.4 10^-3; dt =
3*25.4 10^-3; At = \[Pi]/4 dt^4; Ithroat = (\[Rho]f Lt)/At;
f = 5; n =.; t =.; \[Chi] = 1/f; \[Omega] = 2 \[Pi] f; r = 1.; L =
5 r; \[Theta] = \[Omega] t;
\[Theta]2 = \[Theta] + 2 \[Pi]/3.;
\[Theta]3 = \[Theta]2 + 2 \[Pi]/3;
Clear[y1, y2, s1, P1, P2, P3];
P1[\[Theta]_] := -r Sin[\[Theta]] - (r^2 Sin[\[Theta]] Cos[\[Theta]])/
Sqrt[L^2 - r^2 (Sin[\[Theta]])^2];
P2[\[Theta]_] := P1[\[Theta] + 2 Pi/3];
P3[\[Theta]_] := P1[\[Theta] + 4 Pi/3];
fsum[\[Theta]_] :=
Max[P1[\[Theta]], 0] + Max[P2[\[Theta]], 0] + Max[P3[\[Theta]], 0];
ClearAll[Qin];
Qin[t_?NumericQ] := fsum[\[Theta]];
(*----- system eqns------*)
(* ------ Eqn1: y1 = pf ------*)
Eqn1 = y1'[t] + y1[t] Rf/Ithroat + y2[t] 1 /Ca - Rf Qin[t];
(* ------ Eqn 2: y2 = qf ------*)
Eqn2 = y2'[t] - y1[t] (1/Ithroat);
s1 = NDSolve[{Eqn1 == 0, Eqn2 == 0, y1[0] == 10 10^6,
y2[0.1] == 0.00}, {y1, y2}, {t, 0., tp}, MaxStepSize -> 1/4*tp];
Plot[{Qin[t], Qin[t] - 1/Ithroat Evaluate[y1[t] /. s1]}, {t, 0, 1},
PlotRange -> {{0, 1}, {0, 1.2}}, Frame -> True, Axes -> True,
FrameLabel -> {"Time (secs)", "Output flow rates (L/s)"},
GridLines -> Automatic, AspectRatio -> .65,
PlotStyle -> {{Red, Thickness[0.01]}, {Blue, Thickness[0.01]}},
BaseStyle -> {FontWeight -> "Bold", FontSize -> 15}]