2 Added alternative edited Apr 28 '16 at 13:47 Michael E2 159k1313 gold badges218218 silver badges517517 bronze badges Below is a workaround for the simple case. The OP can say whether it works more general. I haven't quite tracked down yet why the system is set up incorrectly with the default Method -> {"EquationSimplification" -> "Solve"} and with Method -> {"EquationSimplification" -> "Residual"}. But it works in this case with Method -> {"EquationSimplification" -> "MassMatrix"}. simDAE = First@ NDSolve[{x == 100., x'[t] == 0.05 x[t] + deposit[t], y == 1., y[t] == 1. }, {x, y}, {t, 0, 10}, Method -> {"EquationSimplification" -> "MassMatrix"}]; Plot[Evaluate@(x[t] /. simDAE), {t, 0, 10}] Update: Response to comment Again, I can only present a potential workaround at this point. Constructing a WhenEvent[] seems to work better than the automatic processing of DiracDelta[] in this case. For what it's worth, here's a function to convert deposit to a sequence of events, but it's probably easier to use its last line to construct the sequence directly from the times and amounts. ClearAll[diracToEvent]; diracToEvent[depositFn_, x_, t_, scale_: 1/2] := Module[{ times = Union @@ Cases[depositFn, DiracDelta[e_] :> (t /. Solve[e == 0, t]), Infinity], dt, amounts}, With[{xt = If[MatchQ[x, _[t]], x, x[t]]}, dt = Min@Differences@times; amounts = Integrate[depositFn, {t, # - dt*scale, # + dt*scale}] & /@ times; MapThread[ WhenEvent[t > #1, xt -> xt + #2] &, {times, amounts} ] ] ]; simDAE = First@NDSolve[{ x == 100., x'[t] == 0.05 x[t], y[t] == 25 Log[x[t]] , diracToEvent[deposit[t], x, t]}, {x, y}, {t, 0, 10} ]; Plot[Evaluate@({x[t], y[t]} /. simDAE), {t, 0, 10}] Below is a workaround for the simple case. The OP can say whether it works more general. I haven't quite tracked down yet why the system is set up incorrectly with the default Method -> {"EquationSimplification" -> "Solve"} and with Method -> {"EquationSimplification" -> "Residual"}. But it works in this case with Method -> {"EquationSimplification" -> "MassMatrix"}. simDAE = First@ NDSolve[{x == 100., x'[t] == 0.05 x[t] + deposit[t], y == 1., y[t] == 1. }, {x, y}, {t, 0, 10}, Method -> {"EquationSimplification" -> "MassMatrix"}]; Plot[Evaluate@(x[t] /. simDAE), {t, 0, 10}] Below is a workaround for the simple case. The OP can say whether it works more general. I haven't quite tracked down yet why the system is set up incorrectly with the default Method -> {"EquationSimplification" -> "Solve"} and with Method -> {"EquationSimplification" -> "Residual"}. But it works in this case with Method -> {"EquationSimplification" -> "MassMatrix"}. simDAE = First@ NDSolve[{x == 100., x'[t] == 0.05 x[t] + deposit[t], y == 1., y[t] == 1. }, {x, y}, {t, 0, 10}, Method -> {"EquationSimplification" -> "MassMatrix"}]; Plot[Evaluate@(x[t] /. simDAE), {t, 0, 10}] Update: Response to comment Again, I can only present a potential workaround at this point. Constructing a WhenEvent[] seems to work better than the automatic processing of DiracDelta[] in this case. For what it's worth, here's a function to convert deposit to a sequence of events, but it's probably easier to use its last line to construct the sequence directly from the times and amounts. ClearAll[diracToEvent]; diracToEvent[depositFn_, x_, t_, scale_: 1/2] := Module[{ times = Union @@ Cases[depositFn, DiracDelta[e_] :> (t /. Solve[e == 0, t]), Infinity], dt, amounts}, With[{xt = If[MatchQ[x, _[t]], x, x[t]]}, dt = Min@Differences@times; amounts = Integrate[depositFn, {t, # - dt*scale, # + dt*scale}] & /@ times; MapThread[ WhenEvent[t > #1, xt -> xt + #2] &, {times, amounts} ] ] ]; simDAE = First@NDSolve[{ x == 100., x'[t] == 0.05 x[t], y[t] == 25 Log[x[t]] , diracToEvent[deposit[t], x, t]}, {x, y}, {t, 0, 10} ]; Plot[Evaluate@({x[t], y[t]} /. simDAE), {t, 0, 10}] 1 answered Apr 28 '16 at 12:06 Michael E2 159k1313 gold badges218218 silver badges517517 bronze badges Below is a workaround for the simple case. The OP can say whether it works more general. I haven't quite tracked down yet why the system is set up incorrectly with the default Method -> {"EquationSimplification" -> "Solve"} and with Method -> {"EquationSimplification" -> "Residual"}. But it works in this case with Method -> {"EquationSimplification" -> "MassMatrix"}. simDAE = First@ NDSolve[{x == 100., x'[t] == 0.05 x[t] + deposit[t], y == 1., y[t] == 1. }, {x, y}, {t, 0, 10}, Method -> {"EquationSimplification" -> "MassMatrix"}]; Plot[Evaluate@(x[t] /. simDAE), {t, 0, 10}] 