Off[
        General::"spell1", $IterationLimit::"itlim"];
<< Algebra`ReIm`
zrealni := (# /: Im[#] = 0) &;
zrealni /@ {casvypnuti, t};

Vypnutí chybových hlášek, méně zkušeným uživatelům nedoporučujeme.

Dále jsme načetli package, který usnadňuje práci s komplexními čísly a připravili si funkci, která označí svůj argument jako reálné číslo. Mohli jsme ji zavolat na proměnnou casvypnuti; dopředu jsme ale nevěděli, jaké všechny proměnné budeme chtít považovat za reálná čísla a tak ji mapujeme na list, do kterého při psaní programu postupně připisujeme další proměnné.

Rz = (10^3);
Lz = 0.01;
Cz = 10*(10^(-12));
Rt = 5*(10^6);
Lt = 35;
Ct = 15*(10^(-12));
ω = 100*Pi;
T = ((2*Pi)/ω);
Ampl = 10*(10^3);
fi = 0;
Uz = Ampl*Sin[ω*t + fi];

Zadali jsme hodnoty proměnných, které považujeme za známé. Jelikož zadáváme jako hodnoty čísla s nulovou imaginární částí, není nutno na jména těchto proměnných volat funkci "zrej".

par[
          imp_List] := ((((imp^(-1)) /. List -> Plus))^(-1));

Ve schématu se vyskytuje řada paralelních kombinací impedancí. Proto je výhodné zavěst funkci, která listu impedancí přiřadí hodnotu paralení kombinace. Je-li imp={z1,z2,z3}, dá (imp^(-1)) ={(1/z1), (1/z2), (1/z3)}.Tato stuktura vlastně znamená List[(1/z1), (1/z2), (1/z3)].  Chceme-li sečíst např. a+b+c, Je to ekvivalentní struktuře Plus[a,b,c]. Zaměníme-li hlavičku "List"  za hlavičku "Plus", získáme součet, tedy Plus[(1/z1), (1/z2), (1/z3)]=(1/z1) + (1/z2) + (1/z3). Toto umocněno na minus první dá výslednou hodnotu impedance paralelního spojení. Výhodou této syntaxe je její úspornost. Shodný výsledek dá také:
par[imp_List] := Module[{nimp, i}, 
          nimp = Length[imp]; 
          Sum[imp[[i]]^(-1), {i, nimp}]^(-1)].

z1 = 
        par[{I*ω*Lz, Rz}];

    z2 = 
        par[{(1/(I*ω*Cz)), (1/(I*ω*Ct)), I*ω*Lt, 
            Rt}];

    Uf = Ampl*Exp[I*fi];

    Utf = Uf*(z2/(z1 + z2));

    ILtf = (Utf/(I*ω*Lt));

    ILzf = ((Uf - Utf)/(I*ω*Lz));

    vcase[fazor_, cas_] := 
        Im[fazor*Exp[I*ω*cas]];

    iLz0 = vcase[ILzf, casvypnuti];

    iLt0 = vcase[ILtf, casvypnuti];

    uCt0 = vcase[Utf, casvypnuti];

    uCz0 = vcase[Utf, casvypnuti];

K vypnutí dojde -jak předpokládáme- z harmonického ustáleného stavu. Jelikož jde o lineární obvod, řešení najdeme s výhodou pomocí fázorové symboliky. Transformace mezi fázory a harmonickými časovými průběhy je definovába následujícími vztahy:

je-li u[t]=Um*Sin[ω*t+φ], je ospovídající fázor roven Uf=Um*Exp[I*φ]. Fázoru Uf pak odpovídá časový průběh u[t]=Im[Uf*Exp[I*ω*t]]. Pro převod z fázorů do časově závislých veličin jsme si vytvořili funkci "vcase", která má dva parametry: fázor a čas. Pomocí dříve zavedené funkce "par" a vztahů pro dělič napětí získáme fázor napětí na vypínači, když je sepnut. Pro určení počátečních podmínek pro řešení přechodného děje -vypnutí vypínače- použijeme funkci "vcase", ve které za paramery dosadíme příslušný fázor a za čas zatím proměnnou se zatím neurčenou hodnotou "casvypnuti".

rcezdroj = {Uz == Lz*iLz'[t] + uCz[t], 
          iLz[t] + ((Uz - uCz[t])/Rz) == 
            Cz*uCz'[t]};
poczdroj = {iLz[0] == iLz0, 
          uCz[0] == uCz0};
neznamezdroj = {iLz[t], uCz[t]};
reszdroj = DSolve[Union[rcezdroj, poczdroj], neznamezdroj, 
                t][[1]] // Chop // 
          Simplify;
napetizdroj[casvypnuti_] = uCz[t] /. reszdroj;

Po vypnytí se nám obvod rozpadá na dva -až na společnou zem- nezávislé obvody, které popíšeme rovnicemi. Jelikož řešení svěříme Mathematice, není jejich počet rozhodující. Nebudeme tedy psát integro-diferenciální rovnice, ale jakékoli obvodové rovnice, které umíme napsat. Podstatné je, že pro napětí na cívce využijeme vztah uL[t]=L*iL'[t] a pro proud kondenzátorem iC[t]=C*uC'[t]. Velčina, jejíž derivace se v popisu obvodu vyskytne, musí mít určenou také počáteční podmínku. Pro proud ideálním odporem platí rovnice algebraická, tedy uR[t]=R*iR[t], taková rovnice tedy není diferenciální, nezvyšuje řád systému a nezakládá potřebu další počáteční podmínky. Dobré pravidlo je, že volíme takový popis obvodu, kde s uvedenými vztahy vystačíme a využijeme toho, že vždy můžeme popsat obvod algebraickými (bez derivací)  rovnicemi a diferenciálními rovnicemi prvního řádu. Uvedený popis není jediný možný.
Jelikož jsme kromě času sepnutí zadali hodnoty všech proměnných hodnotou a Mathematica tudíž počítá numericky, vycházejí čísla typu 10.030001556505217` - 
        6.447875110334765`*^-22 i. Je zřejmé, že nenulová imaginární část je dána numerickou chybou. Vynulovat taková "skoronulová vůči ostatním" čísla pomůž příkaz Chop.

rcetrafo = {uCt[t] == Lt*iLt'[t], 
          uCt[t]/Rt + Ct*uCt'[t] + iLt[t] == 
            0};
poctrafo = {uCt[0] == uCt0, 
          iLt[0] == iLt0};
neznametrafo = {iLt[t], uCt[t]};
restrafo = DSolve[Union[rcetrafo, poctrafo], neznametrafo, 
                  t][[1]] // Chop // Simplify // 
          N;
napetitrafo[casvypnuti_] = uCt[t] /. restrafo;

Itf = (Utf/z2);
it[t_] = vcase[Itf, t];
tvypopt = 
        ti0 /. NMinimize[{it[ti0]^2, 
                0 ≤ ti0 < 
                  T}, {ti0}][[2]];

Vypnutí znamená, že proud je nulový. Nejpřirozenějším okamžikem je tedy pro vypnutí chvíle, kdy proud nulový je. (Omluvte na tomto místě intuitivní výklad; zamyslete se nad tímto faktem a otázky rádi zodpovíme.) Logické je tedy hledat čas ti0, kdy je proud vypínačem nulový. Znamená to vlastně řešit rovnici: it[ti0]=-0.909119 Cos[100 π ti0]+0.00199967 Sin[100 π ti0]==0.
Nejde o rovnici, která by byla bez použití transcendentních funkcí řešitelná pomocí konečného algoritmu. V Mathematice je možno použít příkaz FindRoot, u kterého je ovšem obtížné řešení omezit do zvoleného intervalu; funkce Sin a Cos jsou ovšem periodické a řešení dané rovnice je tedy nekonečně mnoho. Pro reálné funkce ovšem platí, ře f[xhledane]==0 <=> \(TraditionalForm\`\((f[xhledane])\)\^2\ má\ nulové\ \(\(minimum\)\(.\)\)\) Takže s výhodou použijeme příkaz NMinimize na kvadrát zmíněného proudu; s výhodou, neboť příkaz NMminimize umožňuje zadat boolovskou podmínku, která určí, kde má minimum být. Najdeme optimální čas vypnutí a vložíme jej do proměnné tvypopt.

tmax = T/20;
cvyp = tvypopt;
pts = {RGBColor[0, 0, 1], 
          Thickness[0.006]};
Plot[napetitrafo[cvyp], {t, 0, tmax}, PlotRange -> All, 
        AxesLabel -> {"t[s]", "uT[V]"}, 
        PlotStyle -> pts];
Plot[napetizdroj[cvyp], {t, 0, T}, PlotRange -> All, 
        AxesLabel -> {"t[s]", "uz[V]"}, 
        PlotStyle -> pts];
Plot[napetizdroj[cvyp] - napetitrafo[cvyp], {t, 0, 0.75  T}, 
        PlotRange -> All, AxesLabel -> {"t[s]", "uvyp[V]"}, 
        PlotPoints -> 1000, PlotStyle -> pts];