@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)}] 

@Markiyan Hirnyk 

@Markiyan Hirnyk 

At t= 0  the model is index 1 DAE. That is if you differentiate the second equation you get

2ydy/dt+2zdz/dt=0 

substituting first equation one gets

2zy+2zdz/dt=0

which implies z = 0 or dz/dt = -y

when z is not equal to zero one can solve 2 ODEs

dy/dt = 0 and dz/dt = - y

But when z = 0 we have

dy/dt = 0 and 

y^2- 1 = 0

which has two solutions for y, y =1 and y =-1. We need to differentiate the second equation twice to get an equation for dz/dt (hence this is called index 2 DAE).

When this DAE is sovled in standard index 1 DAE solvers in FORTRAN or MATLAB, the code stops at t = Pi/2 accurately. Maple's dsolve numeric will integrate past this singular point. Code attached for Maple's dsolve numeric. (Uncomment infolevel statement to see the equations being integrated by dsolve/numeric).

Download index2example.mws

@Markiyan Hirnyk 

Markiyan, I use my middle initial as well in publications. My google schoalr page is https://scholar.google.com/citations?user=_rQBjg4AAAAJ&hl=en. The author you are refering to is some one else.

Thanks

Note that when dsolve starts, the time till t = tinit is for initialization. Setting initstep=0.1 (to avoid Maple spending too much time near t =0) and compile =true will help solve many problems efficiently (can't be done for Digits greater than 15 as of Maple 2015).

Modified version of example 4 is attached to show the improvement in the efficiency.

Download example4improved.mws

@acer 

works like a charm in Maple2015. Rosenbrock methods keep the same coefficient matrix at different stages, and this is extremely useful for that.

@roman_pearce 

Roman, I will get back to you on this soon after I gather references. Meanwhile Maple 14 is able to do sparse Linear solve for Digits: = 20 (software float), but Maple 2015 cannot do Linear Solve for more than 15 Digits using the sparse option. Is this because Maple moved to a different sparse solver? It will be great to continue providing high precision solvers.

Thanks

@ 

According to this post from roman_pearce it can't be done. Perhaps Maple will consider providing this option in the future releases