@Axel Vogt 

I think sol2 is for the original ODE and needs the IC. I am using 0.153...as used by the OP.

@Axel Vogt 

Bogo, also for this specific problem, I am not sure if the constants and ICs should be looked at. Both Y(x) and YFD(x) are very close (almost same).

@bogo 

Bogo, you are welcome. Thanks for the acknowledgement. You can call sol1(300) instead of sol1(100) to get sol2 for longer times. Also, for some reason, in Maple 14 classic worksheet it won't plot from 1 to 49 in the same graph, but the code wil plot if you split and plot. You can save those plots separately and display together.

See code below.

Download knownfunctiondsolve3.mws

@bogo 

You are welcome.

I think it is just the syntax to enable dsolve/numeric codes to read. Not sure, perhpas Allan Wittkopf or other Maplesoft developers can explain the need for the same.

I have posted a more efficient approach. See below.

@Preben Alsholm 

Preben, that is perfect. The parameters epsilon(=mu), tinit (initialization time), the constant q inside the tanh function are tunable parameters. My experience suggests that, 

epsilon = 1e-4 is good for abserr =1 e-6,

tinit = 0.1 is a good start for total time of integration of 1. Typically, we don't need more than tinit = 200 if we have to integrate till t = 10,000 for the original DAE system.

initstep = 0.1 or 0.01 is a good number (greater than epsilon).

q= 1000 is good, but for very long integrations, q = 100 may be sufficient.

@Markiyan Hirnyk 

you wrote "You wrote "if you take a linear system of DAEs. The mu approach will theoretically go to the expected steady state irrespective of the value of mu. (When only algebraic equations are solved). For example take Az=b as the algebraic equation". Can you kindly explain that in detail?"

@Markiyan Hirnyk 

not, the modified equations. When I find time, I will post codes for non-autonomous systems that work. If there are specific equations you want me to solve, let me know.

@Markiyan Hirnyk 

df/dt=f(u,z)

0=g(u,z)

rhs of ODEs for non-autonomous system will have f(u,z,t)

@Markiyan Hirnyk 

Sorry for the terminology confusion. Epsilon is used in the pdf file and mu was used in the Maple file.

@Preben Alsholm 

As of now the model written does not have an explicit t dependence in it. Of course one can add a dummy variable dTdt = 1 and convert non-autonomous system to autonomous system. My interests are in large scale systems. So adding one more equation is not bad.

The theory of Rosenbrock methods is well developed, but the implementation in many sovlers do not address this (non-autonomous systems) always.

The proposed approach can handle non-autonomous systems directly by freezing explicit 't' dependent functions in g by multiplying by tanh function. In addition, the explicit time dependent functions in the right hand side (r.h.s.) of f also needs to be taken care of.

We have done this in at least 2-3 different appraoches and are currently trying to find and document the best possible approach by trying different test cases. If you have oscillatory forcing functions, then including partial derivative of that with t may be important.

@Preben Alsholm 

Most of the statements about dsolve are based on standard dsolve approaches inbuilt into Maple. This implicit ODE was given as an example as one can use the proposed approach wihtout using fsolve and including compile=true.

For complicated implicit ODEs, fsolve approach may not work (Because of the insistence on strict tolerance of 1e-15).

Also, please note that first 3 examples can be solved with standard Maple dsolve/numeric command if we use fsolve to find consistent IC for algebraic variable before calling dsolve/numeric. For smaller and linear systems, Maple can extend the Algebraic equations to arrive at the explicit form dz/dt = rhs. However, for nonlinear and large scale DAEs, Maple wastes time and RAM unsuccessfully.

@ 

Acer, is there a reason why we need to enter 0 for zero entries? Is it possible to avoid entering zeros for sparse matrices?

Also, the output  for Solve(3,VV) returns a matrix with 1 column. Is it possible to return a vector?

Thanks

@Markiyan Hirnyk 

if you take a linear system of DAEs. The mu approach will theoretically go to the expected steady state irrespective of the value of mu. (When only algebraic equations are solved). For example take Az=b as the algebraic equation.

@Markiyan Hirnyk 

Markiyan, I modified the text to remove the word bifurcation. The theory of bifurcation for DAEs is more difficult than for ODEs.

I wondered where I saw the word bifurcation. I found the reference. http://link.springer.com/article/10.1007%2FBF01398915

Please check page 555. The authors claim, 

In the first step a simple DAE like (4.1) is analyzed.
(4.1) y'=z, y(0)= 1
y^2 + z^2-1, z(0)=0
Under the index-1 assumption a unique solution exists for t < Pi/2:
(4.2) y (t) = sin (t), z (t) = cos (t)
At t=Pi/2 the index-1 assumption is violated and
(4.3) y(t)= 1, z(t)=0
or bifurcations from (4.3) to (4.2) are admissible solutions.

@Markiyan Hirnyk 

Markiyan,

When I use Maple to solve this I get (Without specifying ICs)

[{y(t) = 1}, {z(t) = 0}], [{y(t) = -1}, {z(t) = 0}], [{y(t) = -sin(-t+_C1), y(t) = sin(-t+_C1)}, {z(t) = diff(y(t),t)}]