Wien oscilator in Matlab and Simulink, problems

So I have this little circuit from which I have to calculate u(t). By hand I have calculated using the Laplace transformations that the correct answer is u(t)=3*sin(t/CR)h(t). I assume the answer here given by my code is halved due to the nature od Dirac delta function in MatLab but I do not know what to do about it. I have also run into a similar problem when I created a simulink model, where i aproximated dirac with a repeated sequence. Any suggestion are welcome.

syms v1 v2 v3 c1 c2 uc1 uc2 Duc1 Duc2 r1 r2 a c ig r
jednacine = [c1*Duc1 + v1/r1 + (v1 - v2)/r2 == 0,...
    v3 - a*v1 ==0,...
    (v2-v1)/r2 - ig + c2*Duc2 == 0,...
    uc2 == v2- v3,...
    uc1 == v1];
jed= eliminate(jednacine, [v1, v2, v3]);
res=solve(jed, [Duc1, Duc2]);
syms uc1(t) uc2(t) ig(t)
jedst=subs([diff(uc1)==res.Duc1;diff(uc2)==res.Duc2],[uc1, uc2, r1, r2, c1, c2, a, ig], [uc1(t), uc2(t), r, r, c, c,3, 1e-6*dirac(t)]);
assume (r>0 & c>0)
resenje=dsolve(jedst, [uc1(0)==0, uc2(0)==0], 'IgnoreAnalyticConstraints', false);
u(t)=subs(3*resenje.uc1, [r,c], [1e3, 1e-6]);
u1(t)=rewrite(simplify(u(t)), 'sin')

u1(t) = 

fplot(t, u1(t), [0, 2e-2])
grid on
xlabel('t [s] ')

[VBBV]

I guess one reason could be when values of both resistor and capacitor (i.e. r1 = r2, c1 = c2) become/are equal, that results in a different expression exactly by half.

Answer 1

[VBBV]

Since you know it's due to nature of Dirac delta function, it would make sense if you modify this line as below

jedst=subs([diff(uc1)==res.Duc1;diff(uc2)==res.Duc2],[uc1, uc2, r1, r2, c1, c2, a, ig], [uc1(t), uc2(t), r, r, c, c,3, 1e-6*2*dirac(t)]);

[Andrea]

Well, it is just an assumption from my side, it is still left unknown why results dont match up.

[VBBV]

syms v1 v2 v3 c1 c2 uc1 uc2 Duc1 Duc2 r1 r2 a c ig r
jednacine = [c1*Duc1 + v1/r1 + (v1 - v2)/r2 == 0,...
    v3 - a*v1 ==0,...
    (v2-v1)/r2 - ig + c2*Duc2 == 0,...
    uc2 == v2- v3,...
    uc1 == v1];
jed= eliminate(jednacine, [v1, v2, v3]);
res=solve(jed, [Duc1, Duc2]);
syms uc1(t) uc2(t) ig(t)
jedst=subs([diff(uc1)==res.Duc1;diff(uc2)==res.Duc2],[uc1, uc2, r1, r2, c1, c2, a, ig], [uc1(t), uc2(t), r, r, c, c,3, 1e-6*2*dirac(t)]);
assume (r>0 & c>0)
resenje=dsolve(jedst, [uc1(0)==0, uc2(0)==0], 'IgnoreAnalyticConstraints', false);
u(t)=subs(3*resenje.uc1, [r,c], [1e3, 1e-6]);
u1(t)=vpa(rewrite(simplify(u(t)), 'sin'),3)

u1(t) = 

fplot(t, u1(t), [0, 2e-2])
grid on
xlabel('t [s] ')