How to compute DSolve or NDSolve for kind of piece wise function?

Question

I need to Compute the expression for two functions. q[t] and q1[t] over certain range of time. The code is given below.

Remove["Global`*"];
e = 1.6*10^-19;
h = 6.63*10^-34;
phi = 0.33*1.6*10^-19;
s = 66*10^-9;
m = 9.11*10^-31;
A = 10^-4;
pi = 3.14;
\[Epsilon]0 = 8.85*10^-12;
k1 = 7;
k2 = 7;
dt = 100*10^-9;
r = 1;
fre = 10000;
v0 = 8.7;
T = .5*(1/fre);
v[t] = (v0*Sin[2*\[Pi]*fre*t]);
ListPlot[Table[{t, %}, {t, 0, T, T/(1000 + 1)}]]

c0 = (A*\[Epsilon]0*k1)/dt;
cap = c0/(1 + (((s*q1[t])/(dt*q[t]))));

eQ = (q[t] + q1[t])/(2*A*\[Epsilon]0*k2);
vk = eQ*s;

ik1 = (e*A)/(2*\[Pi]*h*s^2) ((phi - (e*vk)/2)*E^((-((4*\[Pi]*(Sqrt[2*m])*s)/h))*(Sqrt[phi - (e*vk)/2])) - (phi + (e*vk)/2)*E^((-((4*\[Pi]*(Sqrt[2*m])*s)/h))*(Sqrt[phi + (e*vk)/2])));          (* for evk<phi*)

ik2 = (e^3*A*vk^2)/(4*\[Pi]*h*s^2) (E^(-((4*\[Pi]*(Sqrt[m])*s*phi^1.5)/(h*e*vk))) - (1 + (2*e*vk)/phi)*E^((-((4*\[Pi]*(Sqrt[m])*s*phi^1.5)/(h*e*vk)))*Sqrt[(1 + (2*e*vk)/phi)]));   (* for evk>phi*)                     

(* Below two lines are kind of Pseudocode lines to define "ik" in the dsolve expression . I am so sorry I don't know how to define this kind of if condition in mathematica and hen using in dsolve *) 
ik = If[evk < phi, ik1, ik2];

qq = First@NDSolve[{q'[t] == v[t]/r - q[t]/r*(1/cap), q1'[t] == ik, q[0] == 10^-6, q1[0] == 10^-6}, {q, q1}, t];

Rationalize@ComplexExpand[q[t] /. qq];
Plot[%, {t, 0, T}]
Rationalize@ComplexExpand[q1[t] /. qq];
Plot[%, {t, 0, T}]

I want to calculate the expression for q[t] and q1[t]. The initial conditions can be changed to any value within the range of 10^-12 to 0.

Edit : The Dsolve is not essential, I NDsolve can also be used.

The time range over which NDSolve is DSOLVE is computed can also be change to any range of time.

Answer 1

Results can be obtained with the following changes.

First, delete the two asterisks at the end of the line ik=ik2 (* if evk>phi*)**, which confuse Mathematica.

Second replace DSolve by First@NDSolve (and t by { t, 0, T} in NDSolve), because DSolve cannot solve these equations.

Third, change T to .5*(1/fre) to avoid stiffness issues.

Then, labeling the two curves for convenience,

Corrected Addendum

An earlier version of this addendum produced noisy results for q1[t]. This can be corrected by replacing the various constants in the Question by their exact counterparts

e = 16*10^-20; h = 663*10^-36; phi = (1/3)*16*10^-20; s = 66*10^-9;
m = 911*10^-33; A = 10^-4; ϵ0 = 885*10^-14; k1 = 7; k2 = 7; dt = 100*10^-9;
r = 1; fre = 10000; v0 = 87/10; T = 5*(1/fre);

and increasing WorkingPrecision in NDSolve.

qq = First@NDSolve[{q'[t] == v[t]/r - q[t]/r*(1/cap), q1'[t] == ik, 
       q[0] == 10^-6, q1[0] == 10^-6}, {q, q1}, { t, 0, T}, WorkingPrecision -> 30];

The condition

ik = Piecewise[{{Re[ik1], -phi < e*vk < phi}, {ik2, phi < e*vk}}, 0];

then can be introduced (based on Question 36839) to address the OP's request in the comment below. Doing so yields for q[t]

and q1[t] is constant at 10^-6.

The upper bound of -phi < e*vk < phi on ik1 and the bound phi < e*vkon ik2 in the definition of ik are as requested by the OP. Since ik1 and ik2 are complex for e*vk smaller than -2 phi, some lower bound must be set on the use of ik1, and I chose to set ik equal to zero for -phi < e*vk for lack of a better choice. In fact, the results are insensitive to the precise value of the lower bound, provided that it is not too close to -2 phi. Curiously, ik1 must be wrapped in Abs to avoid errors.

Comment on the use of Piecewise vs If

With the permission of Michael E2, I am reproducing his very use comments below. In general, it's better to use Piecewise rather than If in constructing equations and functions for solvers (NDSolve, Plot, NIntegrate, etc). The solvers have built-in routines for analyzing the discontinuities present in Piecewise. This should be stressed more, since it is not common in other programming languages. The issue is that If uses HoldRest so that ik1 and ik2 are passed to NDSolve unevaluated inside the If. At some point NDSolve substitutes the initial or computed values for q[t], q1[t] into the equations, probably with ReplaceAll. In that case the symbols ik1, ik2 are not altered and afterwards evaluate to contain non-numeric q[t], q1[t]; hence the NDSolve::ndnum error. Piecewise evaluates its arguments.