correction equation in sci_daskr
[scilab.git] / scilab / modules / differential_equations / sci_gateway / cpp / sci_odedc.cpp
1 /*
2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2011 - DIGITEO - Cedric DELAMARRE
4 *
5 * This file must be used under the terms of the CeCILL.
6 * This source file is licensed as described in the file COPYING, which
7 * you should have received as part of this distribution.  The terms
8 * are also available at
9 * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
10 *
11 */
12 /*--------------------------------------------------------------------------*/
13
14 #include "differential_equations_gw.hxx"
15 #include "function.hxx"
16 #include "double.hxx"
17 #include "string.hxx"
18 #include "list.hxx"
19 #include "callable.hxx"
20 #include "differentialequationfunctions.hxx"
21 #include "runvisitor.hxx"
22 #include "context.hxx"
23 #include "checkodeerror.hxx"
24
25 extern "C"
26 {
27 #include "sci_malloc.h"
28 #include "localization.h"
29 #include "Scierror.h"
30 #include "elem_common.h"
31 #include "configvariable_interface.h"
32 #include "sciprint.h"
33 #include "common_structure.h"
34 #include "scifunctions.h"
35 }
36
37 /*--------------------------------------------------------------------------*/
38 types::Function::ReturnValue sci_odedc(types::typed_list &in, int _iRetCount, types::typed_list &out)
39 {
40     // Methode
41     types::String* pStrType     = NULL;
42     const wchar_t* wcsType            = L"lsoda";
43     int meth                    = 0;
44
45     // y0
46     types::Double* pDblY0       = NULL;
47     double* pdYData             = NULL; // contain y0 following by all args data in list case.
48     int sizeOfpdYData           = 0;
49
50     // Other input args
51     types::Double* pDblT0       = NULL;
52     types::Double* pDblT        = NULL;
53     types::Double* pDblRtol     = NULL;
54     types::Double* pDblAtol     = NULL;
55     types::Double* pDblNg       = NULL;
56     types::Double* pDblW        = NULL;
57     types::Double* pDblIw       = NULL;
58     types::Double* pDblNd       = NULL;
59     types::Double* pDblStdel    = NULL;
60
61     // %ODEOPTIONS
62     types::Double* pDblOdeOptions = NULL;
63
64     // Result
65     types::Double* pDblYOut     = NULL;
66
67     // Indicate if the function is given.
68     bool bFuncF     = false;
69     bool bFuncJac   = false;
70     bool bFuncG     = false;
71
72     int iPos        = 0; // Position in types::typed_list "in"
73
74     int sizeYc      = 0;
75     int sizeOfYSize = 1;
76     int* YSize      = NULL;    // YSize(1) = size of y0,
77     // YSize(n) = size of Args(n) in list case.
78
79     C2F(eh0001).mesflg = 1; // flag to control printing of error messages in lapack routine.
80     // 1 means print, 0 means no printing.
81     C2F(eh0001).lunit = 6;  // 6 = stdout
82
83     int one = 1; // used in C2F(dcopy)
84
85     // For root methode
86     int* jroot = NULL;
87
88     // *** check the minimal number of input args. ***
89     if (in.size() < 4)
90     {
91         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "odedc", 6);
92         return types::Function::Error;
93     }
94
95     // *** Get the methode. ***
96     if (in[0]->isString())
97     {
98         pStrType = in[0]->getAs<types::String>();
99         wcsType = pStrType->get(0);
100         iPos++;
101     }
102
103     if (iPos)
104     {
105         if (wcscmp(wcsType, L"adams") == 0)
106         {
107             meth = 1;
108         }
109         else if (wcscmp(wcsType, L"stiff") == 0)
110         {
111             meth = 2;
112         }
113         else if (wcscmp(wcsType, L"root") == 0)
114         {
115             meth = 3;
116         }
117         else if (wcscmp(wcsType, L"discrete") == 0)
118         {
119             meth = 4;
120         }
121         else if (wcscmp(wcsType, L"rk") == 0)
122         {
123             meth = 5;
124         }
125         else if (wcscmp(wcsType, L"rkf") == 0)
126         {
127             meth = 6;
128         }
129         else if (wcscmp(wcsType, L"fix") == 0)
130         {
131             meth = 7;
132         }
133         else
134         {
135             Scierror(999, _("%s: Wrong value for input argument #%d: It must be one of the following strings: adams, stiff, rk, rkf, fix, root or discrete.\n"), "odedc", 1);
136             return types::Function::Error;
137         }
138     }
139
140     // *** check number of output args according the methode. ***
141     if (meth < 3)
142     {
143         if (_iRetCount != 1 && _iRetCount != 3)
144         {
145             Scierror(78, _("%s: Wrong number of output argument(s): %d or %d expected.\n"), "odedc", 1, 3);
146             return types::Function::Error;
147         }
148     }
149     else if (meth == 3)
150     {
151         if (_iRetCount == 3 || _iRetCount > 4)
152         {
153             Scierror(78, _("%s: Wrong number of output argument(s): %d, %d or %d expected.\n"), "odedc", 1, 2, 4);
154             return types::Function::Error;
155         }
156     }
157     else // meth > 3
158     {
159         if (_iRetCount != 1)
160         {
161             Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "odedc", 1);
162             return types::Function::Error;
163         }
164     }
165
166     // *** check type of input args and get it. ***
167     // y0 => [y0c ; y0d]
168     if (in[iPos]->isDouble())
169     {
170         pDblY0 = in[iPos]->getAs<types::Double>();
171         if (pDblY0->isComplex())
172         {
173             Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "odedc", iPos + 1);
174             return types::Function::Error;
175         }
176
177         if (pDblY0->getCols() != 1)
178         {
179             Scierror(999, _("%s: Wrong size for input argument #%d: A real colunm vector expected (n x 1).\n"), "odedc", iPos + 1);
180             return types::Function::Error;
181         }
182     }
183     else
184     {
185         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix\n"), "odedc", iPos + 1);
186         return types::Function::Error;
187     }
188
189     // nd
190     iPos++;
191     if (in[iPos]->isDouble() == false)
192     {
193         Scierror(999, _("%s: Wrong type for input argument #%d: A real scalar expected.\n"), "odedc", iPos + 1);
194         return types::Function::Error;
195     }
196
197     pDblNd = in[iPos]->getAs<types::Double>();
198
199     if (pDblNd->isComplex() || pDblNd->isScalar() == false)
200     {
201         Scierror(999, _("%s: Wrong type for input argument #%d: A real scalar expected.\n"), "odedc", iPos + 1);
202         return types::Function::Error;
203     }
204
205     if (pDblNd->get(0) > pDblY0->getSize())
206     {
207         Scierror(999, _("%s: Wrong value for input argument #%d: Value must not exceeds dimension of argument %d.\n"), "odedc", iPos + 1, iPos);
208         return types::Function::Error;
209     }
210
211     // stdel => hstep and/or delta with delta=0 as default value
212     iPos++;
213     if (in[iPos]->isDouble() == false)
214     {
215         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "odedc", iPos + 1);
216         return types::Function::Error;
217     }
218
219     pDblStdel = in[iPos]->getAs<types::Double>();
220
221     if (pDblStdel->isComplex())
222     {
223         Scierror(999, _("%s: Wrong type for input argument #%d: A real scalar expected.\n"), "odedc", iPos + 1);
224         return types::Function::Error;
225     }
226
227     if (pDblStdel->getSize() > 2)
228     {
229         Scierror(999, _("%s: Wrong size for input argument #%d: %d or %d values expected.\n"), "odedc", iPos + 1, 1, 2);
230         return types::Function::Error;
231     }
232
233     // t0
234     iPos++;
235     if (in[iPos]->isDouble() == false)
236     {
237         Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "odedc", iPos + 1);
238         return types::Function::Error;
239     }
240
241     pDblT0 = in[iPos]->getAs<types::Double>();
242
243     if (pDblT0->isScalar() == false)
244     {
245         Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "odedc", iPos + 1);
246         return types::Function::Error;
247     }
248
249     // t
250     iPos++;
251     if (in[iPos]->isDouble() == false)
252     {
253         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "odedc", iPos + 1);
254         return types::Function::Error;
255     }
256
257     pDblT = in[iPos]->getAs<types::Double>();
258
259     // size of continuous part
260     sizeYc = pDblY0->getSize() - (int)pDblNd->get(0);
261     if (sizeYc < 0)
262     {
263         Scierror(999, _("%s: Wrong value for input argument #%d: Value of nd exceeds dimension of y0.\n"), "odedc", 2);
264         return types::Function::Error;
265     }
266
267     // get next inputs
268     DifferentialEquationFunctions* deFunctionsManager = new DifferentialEquationFunctions(L"odedc");
269     DifferentialEquation::addDifferentialEquationFunctions(deFunctionsManager);
270     deFunctionsManager->setOdedcYDSize((int)pDblNd->get(0));
271     deFunctionsManager->setOdeYRows(pDblY0->getRows());
272     deFunctionsManager->setOdeYCols(pDblY0->getCols());
273
274     YSize = (int*)MALLOC(sizeOfYSize * sizeof(int));
275     *YSize = pDblY0->getSize();
276     pdYData = (double*)MALLOC(pDblY0->getSize() * sizeof(double));
277     C2F(dcopy)(YSize, pDblY0->get(), &one, pdYData, &one);
278
279     if (meth == 4)
280     {
281         if (in.size() != 7)
282         {
283             Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "odedc", 7);
284             DifferentialEquation::removeDifferentialEquationFunctions();
285             FREE(pdYData);
286             FREE(YSize);
287             return types::Function::Error;
288         }
289
290         if (in[6]->isCallable() == false && in[6]->isString() == false && in[6]->isList() == false)
291         {
292             Scierror(999, _("%s: Wrong type for input argument #%d: A function expected.\n"), "odedc", 7);
293             DifferentialEquation::removeDifferentialEquationFunctions();
294             FREE(pdYData);
295             FREE(YSize);
296             return types::Function::Error;
297         }
298     }
299
300     for (iPos++; iPos < in.size(); iPos++)
301     {
302         if (in[iPos]->isDouble())
303         {
304             if (pDblRtol == NULL && bFuncF == false)
305             {
306                 pDblRtol = in[iPos]->getAs<types::Double>();
307                 if (pDblRtol->getSize() != pDblY0->getSize() && pDblRtol->isScalar() == false)
308                 {
309                     Scierror(267, _("%s: Arg %d and arg %d must have equal dimensions.\n"), "odedc", pStrType ? 2 : 1, iPos + 1);
310                     DifferentialEquation::removeDifferentialEquationFunctions();
311                     FREE(pdYData);
312                     FREE(YSize);
313                     return types::Function::Error;
314                 }
315             }
316             else if (pDblAtol == NULL && bFuncF == false)
317             {
318                 pDblAtol = in[iPos]->getAs<types::Double>();
319                 if (pDblAtol->getSize() != pDblY0->getSize() && pDblAtol->isScalar() == false)
320                 {
321                     Scierror(267, _("%s: Arg %d and arg %d must have equal dimensions.\n"), "odedc", pStrType ? 2 : 1, iPos + 1);
322                     DifferentialEquation::removeDifferentialEquationFunctions();
323                     FREE(pdYData);
324                     FREE(YSize);
325                     return types::Function::Error;
326                 }
327             }
328             else if (pDblNg == NULL && bFuncF == true && meth == 3)
329             {
330                 pDblNg = in[iPos]->getAs<types::Double>();
331             }
332             else if (pDblW == NULL && bFuncF == true && (bFuncG == true || meth != 3))
333             {
334                 if (in.size() == iPos + 2)
335                 {
336                     if (in[iPos + 1]->isDouble() == false)
337                     {
338                         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "odedc", iPos + 2);
339                         DifferentialEquation::removeDifferentialEquationFunctions();
340                         FREE(pdYData);
341                         FREE(YSize);
342                         return types::Function::Error;
343                     }
344
345                     pDblW = in[iPos]->getAs<types::Double>();
346                     pDblIw = in[iPos + 1]->getAs<types::Double>();
347                     iPos++;
348                 }
349                 else
350                 {
351                     Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "odedc", iPos + 2);
352                     DifferentialEquation::removeDifferentialEquationFunctions();
353                     FREE(pdYData);
354                     FREE(YSize);
355                     return types::Function::Error;
356                 }
357             }
358             else
359             {
360                 Scierror(999, _("%s: Wrong type for input argument #%d: A function expected.\n"), "odedc", iPos + 1);
361                 DifferentialEquation::removeDifferentialEquationFunctions();
362                 FREE(pdYData);
363                 FREE(YSize);
364                 return types::Function::Error;
365             }
366         }
367         else if (in[iPos]->isCallable())
368         {
369             types::Callable* pCall = in[iPos]->getAs<types::Callable>();
370             if (bFuncF == false)
371             {
372                 deFunctionsManager->setFFunction(pCall);
373                 bFuncF = true;
374             }
375             else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
376             {
377                 deFunctionsManager->setJacFunction(pCall);
378                 bFuncJac = true;
379             }
380             else if (bFuncG == false && meth == 3)
381             {
382                 deFunctionsManager->setGFunction(pCall);
383                 bFuncG = true;
384             }
385             else
386             {
387                 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "odedc", iPos + 1);
388                 DifferentialEquation::removeDifferentialEquationFunctions();
389                 FREE(pdYData);
390                 FREE(YSize);
391                 return types::Function::Error;
392             }
393         }
394         else if (in[iPos]->isString())
395         {
396             types::String* pStr = in[iPos]->getAs<types::String>();
397             bool bOK = false;
398
399             if (bFuncF == false)
400             {
401                 bOK = deFunctionsManager->setFFunction(pStr);
402                 bFuncF = true;
403             }
404             else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
405             {
406                 bOK = deFunctionsManager->setJacFunction(pStr);
407                 bFuncJac = true;
408             }
409             else if (bFuncG == false && meth == 3)
410             {
411                 bOK = deFunctionsManager->setGFunction(pStr);
412                 bFuncG = true;
413             }
414             else
415             {
416                 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "odedc", iPos + 1);
417                 DifferentialEquation::removeDifferentialEquationFunctions();
418                 FREE(pdYData);
419                 FREE(YSize);
420                 return types::Function::Error;
421             }
422
423             if (bOK == false)
424             {
425                 char* pst = wide_string_to_UTF8(pStr->get(0));
426                 Scierror(50, _("%s: Subroutine not found: %s\n"), "odedc", pst);
427                 FREE(pst);
428                 DifferentialEquation::removeDifferentialEquationFunctions();
429                 FREE(pdYData);
430                 FREE(YSize);
431                 return types::Function::Error;
432             }
433         }
434         else if (in[iPos]->isList())
435         {
436             types::List* pList = in[iPos]->getAs<types::List>();
437
438             if (pList->getSize() == 0)
439             {
440                 Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "odedc", iPos + 1, "(string empty)");
441                 DifferentialEquation::removeDifferentialEquationFunctions();
442                 FREE(pdYData);
443                 FREE(YSize);
444                 return types::Function::Error;
445             }
446
447             if (bFuncF && (bFuncJac || pDblNg) && (bFuncG || meth != 3))
448             {
449                 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "odedc", iPos + 1);
450                 DifferentialEquation::removeDifferentialEquationFunctions();
451                 FREE(pdYData);
452                 FREE(YSize);
453                 return types::Function::Error;
454             }
455
456             if (pList->get(0)->isString())
457             {
458                 types::String* pStr = pList->get(0)->getAs<types::String>();
459                 bool bOK = false;
460
461                 if (bFuncF == false)
462                 {
463                     bFuncF = true;
464                     bOK = deFunctionsManager->setFFunction(pStr);
465                     sizeOfpdYData = *YSize;
466                 }
467                 else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
468                 {
469                     bFuncJac = true;
470                     bOK = deFunctionsManager->setJacFunction(pStr);
471                     if (sizeOfpdYData == 0)
472                     {
473                         sizeOfpdYData = *YSize;
474                     }
475                 }
476                 else if (bFuncG == false && meth == 3)
477                 {
478                     bFuncG = true;
479                     bOK = deFunctionsManager->setGFunction(pStr);
480                     if (sizeOfpdYData == 0)
481                     {
482                         sizeOfpdYData = *YSize;
483                     }
484                 }
485
486                 if (bOK == false)
487                 {
488                     char* pst = wide_string_to_UTF8(pStr->get(0));
489                     Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "odedc", iPos + 1, pst);
490                     FREE(pst);
491                     DifferentialEquation::removeDifferentialEquationFunctions();
492                     FREE(pdYData);
493                     FREE(YSize);
494                     return types::Function::Error;
495                 }
496
497                 // store in YSize the size of "args" input
498                 int* sizeTemp = YSize;
499                 int totalSize = sizeOfpdYData;
500
501                 YSize = (int*)MALLOC((sizeOfYSize + pList->getSize() - 1) * sizeof(int));
502                 memcpy(YSize, sizeTemp, sizeOfYSize * sizeof(int));
503
504                 std::vector<types::Double*> vpDbl;
505                 for (int iter = 0; iter < pList->getSize() - 1; iter++)
506                 {
507                     if (pList->get(iter + 1)->isDouble() == false)
508                     {
509                         Scierror(999, _("%s: Wrong type for input argument #%d: Argument %d in the list must be a matrix.\n"), "odedc", iPos + 1, iter + 1);
510                         FREE(pdYData);
511                         FREE(YSize);
512                         return types::Function::Error;
513                     }
514
515                     vpDbl.push_back(pList->get(iter + 1)->getAs<types::Double>());
516                     YSize[sizeOfYSize + iter] = vpDbl[iter]->getSize();
517                     totalSize += YSize[sizeOfYSize + iter];
518                 }
519
520                 // store in pdYData "args" input
521                 double* pdYDataTemp = pdYData;
522                 pdYData = (double*)MALLOC(totalSize * sizeof(double));
523                 C2F(dcopy)(&sizeOfpdYData, pdYDataTemp, &one, pdYData, &one);
524
525                 int position = sizeOfpdYData;
526                 for (int iter = 0; iter < pList->getSize() - 1; iter++)
527                 {
528                     C2F(dcopy)(&YSize[sizeOfYSize + iter], vpDbl[iter]->get(), &one, &pdYData[position], &one);
529                     position += vpDbl[iter]->getSize();
530                 }
531                 vpDbl.clear();
532                 sizeOfpdYData = totalSize;
533                 sizeOfYSize += pList->getSize() - 1;
534                 FREE(pdYDataTemp);
535                 FREE(sizeTemp);
536             }
537             else if (pList->get(0)->isCallable())
538             {
539                 if (bFuncF == false)
540                 {
541                     bFuncF = true;
542                     deFunctionsManager->setFFunction(pList->get(0)->getAs<types::Callable>());
543                     for (int iter = 1; iter < pList->getSize(); iter++)
544                     {
545                         deFunctionsManager->setFArgs(pList->get(iter)->getAs<types::InternalType>());
546                     }
547                 }
548                 else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
549                 {
550                     bFuncJac = true;
551                     deFunctionsManager->setJacFunction(pList->get(0)->getAs<types::Callable>());
552                     for (int iter = 1; iter < pList->getSize(); iter++)
553                     {
554                         deFunctionsManager->setJacArgs(pList->get(iter)->getAs<types::InternalType>());
555                     }
556                 }
557                 else if (bFuncG == false && meth == 3)
558                 {
559                     bFuncG = true;
560                     deFunctionsManager->setGFunction(pList->get(0)->getAs<types::Callable>());
561                     for (int iter = 1; iter < pList->getSize(); iter++)
562                     {
563                         deFunctionsManager->setGArgs(pList->get(iter)->getAs<types::InternalType>());
564                     }
565                 }
566             }
567             else
568             {
569                 Scierror(999, _("%s: Wrong type for input argument #%d: The first argument in the list must be a string or a function.\n"), "odedc", iPos + 1);
570                 DifferentialEquation::removeDifferentialEquationFunctions();
571                 FREE(pdYData);
572                 FREE(YSize);
573                 return types::Function::Error;
574             }
575         }
576         else
577         {
578             Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or a function expected.\n"), "odedc", iPos + 1);
579             DifferentialEquation::removeDifferentialEquationFunctions();
580             FREE(pdYData);
581             FREE(YSize);
582             return types::Function::Error;
583         }
584     }
585
586     if (bFuncF == false)
587     {
588         int val = (meth == 3) ? 3 : 1;
589         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "odedc", in.size() + val);
590         DifferentialEquation::removeDifferentialEquationFunctions();
591         FREE(pdYData);
592         FREE(YSize);
593         return types::Function::Error;
594     }
595     if (pDblNg == NULL && meth == 3)
596     {
597         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "odedc", in.size() + 2);
598         DifferentialEquation::removeDifferentialEquationFunctions();
599         FREE(pdYData);
600         FREE(YSize);
601         return types::Function::Error;
602     }
603     if (bFuncG == false && meth == 3)
604     {
605         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "odedc", in.size() + 1);
606         DifferentialEquation::removeDifferentialEquationFunctions();
607         FREE(pdYData);
608         FREE(YSize);
609         return types::Function::Error;
610     }
611
612     // *** Initialization. ***
613
614     int itol        = 1;
615     int itask       = 4;
616     int istate      = 1; // istate = integer flag (input and output).  set istate = 1.
617     int iopt        = 0; // iopt   = 0 to indicate no optional inputs used.
618     int jt          = 2; // jacobian type indicator.  set jt = 2.
619     int ml          = -1;
620     int mu          = -1;
621
622     // work tab
623     double* rwork   = NULL;
624     int* iwork      = NULL;
625     int rworkSize   = 0;
626     int iworkSize   = 0;
627
628     // contain ls0001, lsa001 and eh0001 structures
629     double* dStructTab  = NULL;
630     int* iStructTab     = NULL;
631     int dStructTabSize  = 0;
632     int iStructTabSize  = 0;
633
634     int rwSize  = 0; // rwSize = dStructTab + rworkSize
635     int iwSize  = 0; // iwSize = iStructTab + iworkSize
636
637     // structures used by lsoda and lsode
638     double* ls0001d = &(C2F(ls0001).tret);
639     int* ls0001i    = &(C2F(ls0001).illin);
640     double* lsa001d = &(C2F(lsa001).tsw);
641     int* lsa001i    = &(C2F(lsa001).insufr);
642     double* lsr001d = &(C2F(lsr001).rownr3[0]);
643     int* lsr001i    = &(C2F(lsr001).lg0);
644     int* eh0001i    = &(C2F(eh0001).mesflg);
645
646     // get %ODEOPTIONS
647     types::InternalType* pIT = symbol::Context::getInstance()->get(symbol::Symbol(L"%ODEOPTIONS"));
648     if (pIT != NULL && pIT->isDouble())
649     {
650         pDblOdeOptions = pIT->getAs<types::Double>();
651         if (pDblOdeOptions->getSize() == 12)
652         {
653             iopt    = 1;
654             itask   = (int)pDblOdeOptions->get(0);
655             jt      = (int)pDblOdeOptions->get(5);
656             ml      = (int)pDblOdeOptions->get(10);
657             mu      = (int)pDblOdeOptions->get(11);
658
659             if (itask != 4 && itask != 5)
660             {
661                 if (itask == 1)
662                 {
663                     itask = 4;
664                 }
665                 else
666                 {
667                     itask = 5;
668                 }
669
670                 if (getWarningMode())
671                 {
672                     sciprint(_("%s: Warning: odedc forces itask = %d.\n"), "odedc", itask);
673                 }
674             }
675         }
676     }
677
678     if (iopt == 1 && (pDblOdeOptions->get(4) > pDblOdeOptions->get(3))) // hmin > hmax ?
679     {
680         Scierror(9999, _("%s: Wrong value of hmin and hmax: hmin = %d is greater than hmax = %d.\n"), "odedc", pDblOdeOptions->get(4), pDblOdeOptions->get(3));
681         DifferentialEquation::removeDifferentialEquationFunctions();
682         FREE(pdYData);
683         FREE(YSize);
684         return types::Function::Error;
685     }
686
687     if (jt < 0 || jt > 5)
688     {
689         Scierror(9999, _("%s: Wrong value of Jacobian type: A number between %d and %d expected.\n"), "odedc", 0, 5);
690         DifferentialEquation::removeDifferentialEquationFunctions();
691         FREE(pdYData);
692         FREE(YSize);
693         return types::Function::Error;
694     }
695
696     if (iopt == 0 && bFuncJac)
697     {
698         jt = 1;
699     }
700
701     if (bFuncJac && (jt == 2 || jt == 5) && getWarningMode())
702     {
703         sciprint(_("%s: Warning: Jacobian is given, but not used.\n"), "odedc");
704     }
705
706     if (bFuncJac == false && (jt == 1 || jt == 4))
707     {
708         if (getWarningMode())
709         {
710             sciprint(_("%s: Warning: No Jacobian external given, but one is required by %ODEOPTIONS(6) value. Jacobian will be estimated.\n"), "odedc");
711         }
712
713         jt = 2;
714     }
715
716     //compute itol and set the tolerances rtol and atol.
717     double* rtol = NULL;
718     double* atol = NULL;
719
720     if (pDblRtol)
721     {
722         if (pDblRtol->isScalar())
723         {
724             rtol = (double*)MALLOC(sizeof(double));
725             *rtol = pDblRtol->get(0);
726         }
727         else
728         {
729             rtol = pDblRtol->get();
730             itol += 2;
731         }
732     }
733     else
734     {
735         rtol = (double*)MALLOC(sizeof(double));
736         if (meth == 6 || meth == 7)
737         {
738             *rtol = 1.e-3;
739         }
740         else
741         {
742             *rtol = 1.0e-7;
743         }
744     }
745
746     if (pDblAtol)
747     {
748         if (pDblAtol->isScalar())
749         {
750             atol = (double*)MALLOC(sizeof(double));
751             *atol = pDblAtol->get(0);
752         }
753         else
754         {
755             atol = pDblAtol->get();
756             itol ++;
757         }
758     }
759     else
760     {
761         atol = (double*)MALLOC(sizeof(double));
762         if (meth == 6 || meth == 7)
763         {
764             *atol = 1.e-4;
765         }
766         else
767         {
768             *atol = 1.0e-9;
769         }
770     }
771
772     // Compute rwork, iwork size and create them.
773     int nyh     = (sizeYc);
774     int mxordn  = 12;
775     int mxords  = 5;
776     int lmat    = 0;
777     int lrn     = 0;
778     int lrs     = 0;
779
780     if (pDblW) // structure ls0001 have been restored.
781     {
782         nyh = C2F(ls0001).nyh;
783     }
784
785     if (iopt == 1)
786     {
787         mxordn = (int)pDblOdeOptions->get(7);
788         mxords = (int)pDblOdeOptions->get(8);
789     }
790
791     if (mxordn > 12 || mxords > 5 || mxordn < 1 || mxords < 1)
792     {
793         if (getWarningMode())
794         {
795             sciprint(_("%s: Warning: Wrong value for maximun stiff/non-stiff order allowed :\nAt most %d for mxordn, %d for mxords and no null value for both expected.\nWrong value will be reduced to the default value.\n"), "ode", 12, 5);
796         }
797
798         mxordn = 12;
799         mxords = 5;
800     }
801
802     int maxord = mxords;
803
804     switch (meth)
805     {
806         case 3: // lsodar (root)
807         {
808             //             lrn = 20 + nyh*(mxordn+1) + 3*neq + 3*ng,
809             //             lrs = 20 + nyh*(mxords+1) + 3*neq + lmat + 3*ng,
810             //          where
811             //             nyh    = the initial value of neq,
812             //             mxordn = 12, unless a smaller value is given as an
813             //                      optional input,
814             //             mxords = 5, unless a smaller value is given as an
815             //                      optional input,
816             //             lmat   = length of matrix work space..
817             //             lmat   = neq**2 + 2              if jt = 1 or 2,
818             //             lmat   = (2*ml + mu + 1)*neq + 2 if jt = 4 or 5.
819
820             lrn = 3 * (int)pDblNg->get(0);
821             lrs = lrn;
822             dStructTabSize = 246 - 241;
823             iStructTabSize = 59 - 50;
824             jroot = (int*)MALLOC((int)pDblNg->get(0) * sizeof(int));
825         }
826         case 0: // lsoda
827         {
828             //             lrn = 20 + nyh*(mxordn+1) + 3*neq,
829             //             lrs = 20 + nyh*(mxords+1) + 3*neq + lmat,
830             //          where
831             //             nyh    = the initial value of neq,
832             //             mxordn = 12, unless a smaller value is given as an
833             //                      optional input,
834             //             mxords = 5, unless a smaller value is given as an
835             //                      optional input,
836             //             lmat   = length of matrix work space..
837             //             lmat   = neq**2 + 2              if jt = 1 or 2,
838             //             lmat   = (2*ml + mu + 1)*neq + 2 if jt = 4 or 5.
839
840             if (jt == 1 || jt == 2)
841             {
842                 lmat = (sizeYc) * (sizeYc) + 2;       //  if jt = 1 or 2,
843             }
844             else if (jt == 4 || jt == 5)
845             {
846                 //ml and mu = -1 in all cases
847                 lmat = (2 * ml + mu + 1) * (sizeYc) + 2; //  if jt = 4 or 5.
848             }
849
850             lrn += 20 + nyh * (mxordn + 1) + 3 * (sizeYc);
851             lrs += 20 + nyh * (mxords + 1) + 3 * (sizeYc) + lmat;
852
853             rworkSize   = max(lrn, lrs);
854             iworkSize   = 20 + sizeYc;
855
856             dStructTabSize += 241;
857             iStructTabSize += 50;
858
859             break;
860         }
861         case 1: // lsode (adams) (non stiff)
862         {
863             maxord = mxordn;
864         }
865         case 2: // lsode (stiff)
866         {
867             //          20 + nyh*(maxord + 1) + 3*neq + lmat
868             //        where
869             //          nyh    = the initial value of neq,
870             //          maxord = 12 (if meth = 1) or 5 (if meth = 2) (unless a
871             //                   smaller value is given as an optional input),
872             //          lmat   = 0             if miter = 0,
873             //          lmat   = neq**2 + 2    if miter is 1 or 2,
874             //          lmat   = neq + 2       if miter = 3, and
875             //          lmat   = (2*ml+mu+1)*neq + 2 if miter is 4 or 5.
876
877             if (jt == 1 || jt == 2)
878             {
879                 lmat = (sizeYc) * (sizeYc) + 2;
880             }
881             else if (jt == 3)
882             {
883                 lmat = (sizeYc) + 2;
884             }
885             else if (jt == 4 || jt == 5)
886             {
887                 lmat = (2 * ml + mu + 1) * (sizeYc) + 2;
888             }
889
890             rworkSize = 20 + nyh * (maxord + 1) + 3 * (sizeYc) + lmat;
891             iworkSize = 20; // if jt = 0 or 3
892             if (jt == 1 || jt == 2 || jt == 4 || jt == 5) // iSize = 20 + neq  otherwise
893             {
894                 iworkSize += (sizeYc);
895             }
896
897             dStructTabSize = 219;
898             iStructTabSize = 41;
899             break;
900         }
901         case 4: // lsdisc (discrete)
902         {
903             rworkSize = sizeYc;
904             iworkSize = 1;
905             break;
906         }
907         case 5: // lsrgk (rk)
908         {
909             rworkSize = 3 * (sizeYc);
910             iworkSize = 1;
911             break;
912         }
913         case 6: // rkf45 (rkf)
914         {
915             rworkSize = 3 + 8 * (sizeYc);
916             iworkSize = 5;
917             break;
918         }
919         case 7: // rksimp (fix)
920         {
921             rworkSize = 3 + 8 * (sizeYc);
922             iworkSize = 1;
923             break;
924         }
925     }
926
927     rwSize = rworkSize + dStructTabSize;
928     iwSize = iworkSize + iStructTabSize;
929
930     rwork = (double*)MALLOC(rworkSize * sizeof(double));
931     iwork = (int*)MALLOC(iworkSize * sizeof(int));
932
933     if (meth < 4)
934     {
935         if (pDblW && pDblIw)
936         {
937             if (pDblW->getSize() != rwSize || pDblIw->getSize() != iwSize)
938             {
939                 Scierror(9999, _("%s: Wrong size for w and iw: w = %d and iw = %d expected.\n"), "odedc", rwSize, iwSize);
940                 DifferentialEquation::removeDifferentialEquationFunctions();
941                 FREE(pdYData);
942                 FREE(YSize);
943                 if (jroot)
944                 {
945                     FREE(jroot);
946                 }
947                 if (itol == 1 || itol == 3)
948                 {
949                     FREE(atol);
950                 }
951                 if (itol < 3)
952                 {
953                     FREE(rtol);
954                 }
955                 return types::Function::Error;
956             }
957
958             istate = 2; // 1 means this is the first call | 2  means this is not the first call
959
960             dStructTab = (double*)MALLOC(dStructTabSize * sizeof(double));
961             iStructTab = (int*)MALLOC(iStructTabSize * sizeof(int));
962
963             C2F(dcopy)(&rworkSize, pDblW->get(), &one, rwork, &one);
964             C2F(dcopy)(&dStructTabSize, pDblW->get() + rworkSize, &one, dStructTab, &one);
965
966             for (int i = 0; i < iworkSize; i++)
967             {
968                 iwork[i] = (int)pDblIw->get(i);
969             }
970
971             for (int i = 0; i < iStructTabSize; i++)
972             {
973                 iStructTab[i] = (int)pDblIw->get(i + iworkSize);
974             }
975         }
976         else
977         {
978             //if iopt = 0 default value are used without set it.
979             if (iopt == 1)
980             {
981                 if (itask >= 4)
982                 {
983                     rwork[0] = pDblOdeOptions->get(1); // tcrit
984                 }
985
986                 rwork[4] = pDblOdeOptions->get(2);      // h0 =0
987                 rwork[5] = pDblOdeOptions->get(3);      // hmax = %inf
988                 rwork[6] = pDblOdeOptions->get(4);      // hmin = 0
989                 iwork[0] = (int)pDblOdeOptions->get(10);// ml = -1
990                 iwork[1] = (int)pDblOdeOptions->get(11);// mu = -1
991                 iwork[4] = (int)pDblOdeOptions->get(9); // ixpr = 0
992                 iwork[5] = (int)pDblOdeOptions->get(6); // mxstep = 500
993                 iwork[6] = 10;  // mxhnil = 10 maximum number of messages printed per problem
994                 iwork[7] = (int)pDblOdeOptions->get(7); // mxordn = 12
995                 iwork[8] = (int)pDblOdeOptions->get(8); // mxords = 5
996             }
997         }
998     }
999
1000     if (pDblOdeOptions && pDblOdeOptions->get(9) == 1)
1001     {
1002         sciprint(_("itask = %d\tmeth = %d\tjactyp = %d\tml = %d\tmu = %d\n"), itask, meth, jt, ml, mu, iopt);
1003         sciprint(_("tcrit = %lf\th0 = %lf\thmax = %lf\thmin = %lf\n"), pDblOdeOptions->get(1), pDblOdeOptions->get(2), pDblOdeOptions->get(3), pDblOdeOptions->get(4));
1004     }
1005
1006     // get rwork and iwork
1007     if (meth < 4)
1008     {
1009         if (pDblW && pDblIw)
1010         {
1011             int dPos    = 0;
1012             int iPos    = 0;
1013             int dSize   = 219;
1014
1015             //restore ls0001 and eh0001 from w (rwork) and iw (iwork).
1016             C2F(dcopy)(&dSize, dStructTab, &one, ls0001d, &one);
1017             memcpy(ls0001i, iStructTab, 39 * sizeof(int));
1018
1019             dPos = dSize;
1020             iPos = 39;
1021
1022             //restore lsa001 from w (rwork) and iw (iwork).
1023             if (meth == 0 || meth == 3)
1024             {
1025                 dSize = 22;
1026                 C2F(dcopy)(&dSize, dStructTab + dPos, &one, lsa001d, &one);
1027                 memcpy(lsa001i, &iStructTab[iPos], 9 * sizeof(int));
1028
1029                 dPos += dSize;
1030                 iPos += 9;
1031             }
1032
1033             //restore lsr001
1034             if (meth == 3)
1035             {
1036                 dSize = 5;
1037                 C2F(dcopy)(&dSize, dStructTab + dPos, &one, lsr001d, &one);
1038                 memcpy(lsr001i, &iStructTab[iPos], 9 * sizeof(int));
1039             }
1040         }
1041     }
1042
1043     // *** Perform operation. ***
1044     double nhpass   = 0;
1045     double delta    = 0;
1046     double tleft    = pDblT0->get(0);
1047     double tright   = pDblT0->get(0);
1048     double hf       = 0;
1049     double tf       = 0;
1050
1051     bool bStore     = false;
1052     bool bUpdate    = false;
1053     bool bOneStep   = false;
1054
1055     std::string strMeth;
1056     int err = 0;
1057
1058     if (pDblStdel->isScalar() == false)
1059     {
1060         delta = (int)pDblStdel->get(1);
1061     }
1062
1063     if (itask == 5 || itask == 3 || itask == 2)
1064     {
1065         bOneStep = true;
1066         if (getWarningMode() && pDblT->isScalar() == false)
1067         {
1068             sciprint(_("itask = %d: At most one value of t is allowed, the last element of t is used.\n"), itask);
1069         }
1070     }
1071
1072     if (bOneStep)
1073     {
1074         std::list<double*> pDblYOutList = std::list<double*>();
1075         std::list<double> pDblTOutList = std::list<double>();
1076         int iLastT = pDblT->getSize() - 1;
1077         bool bIntegrateContPart = false;
1078
1079         do
1080         {
1081             if (bIntegrateContPart == false)
1082             {
1083                 hf = min(pDblT0->get(0) + (nhpass + delta) * pDblStdel->get(0), pDblT->get(iLastT));
1084             }
1085
1086             if (fabs(tleft - hf) < 1.0e-12) // update discrete part
1087             {
1088                 bIntegrateContPart = false;
1089                 deFunctionsManager->setOdedcFlag();
1090
1091                 if (pDblOdeOptions && pDblOdeOptions->get(9) == 1)
1092                 {
1093                     sciprint(_("update at t = %lf\n"), tright);
1094                 }
1095                 try
1096                 {
1097                     ode_f(&sizeYc, &tright, pdYData, pdYData + sizeYc);
1098                 }
1099                 catch (ast::ScilabError &e)
1100                 {
1101                     char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
1102                     sciprint(_("%s: Update failed at t = %lf\n"), "odedc", tright);
1103                     Scierror(999, pstrMsg);
1104
1105                     DifferentialEquation::removeDifferentialEquationFunctions();
1106                     FREE(pdYData);
1107                     FREE(YSize);
1108                     if (jroot)
1109                     {
1110                         FREE(jroot);
1111                     }
1112                     if (dStructTab)
1113                     {
1114                         FREE(dStructTab);
1115                     }
1116                     if (iStructTab)
1117                     {
1118                         FREE(iStructTab);
1119                     }
1120                     if (itol == 1 || itol == 3)
1121                     {
1122                         FREE(atol);
1123                     }
1124                     if (itol < 3)
1125                     {
1126                         FREE(rtol);
1127                     }
1128                     return types::Function::Error;
1129                 }
1130
1131                 deFunctionsManager->resetOdedcFlag();
1132                 nhpass++;
1133
1134                 double* copy = (double*)MALLOC(*YSize * sizeof(double));
1135                 C2F(dcopy)(YSize, pdYData, &one, copy, &one);
1136                 pDblYOutList.push_back(copy);
1137                 pDblTOutList.push_back(tleft);
1138             }
1139             else // integrate continuous part
1140             {
1141                 bIntegrateContPart = true;
1142                 tright = hf;
1143                 rwork[0] = hf; // tcrit = hf
1144
1145                 if (pDblOdeOptions && pDblOdeOptions->get(9) == 1)
1146                 {
1147                     sciprint(_("integ. from tleft= %lf to tf= %lf\n"), tleft, tright);
1148                 }
1149
1150                 try
1151                 {
1152                     switch (meth)
1153                     {
1154                         case 1: // lsode (adams)
1155                         case 2: // lsode (stiff)
1156                         {
1157                             strMeth = "lsode";
1158                             int jacType = 10 * meth + jt;
1159                             C2F(lsode)(ode_f, &sizeYc, pdYData, &tleft, &tright, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jacType);
1160                             break;
1161                         }
1162                         case 3: // lsodar
1163                         {
1164                             strMeth = "lsodar";
1165                             int ng = (int)pDblNg->get(0);
1166                             C2F(lsodar)(ode_f, &sizeYc, pdYData, &tleft, &tright, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jt, ode_g, &ng, jroot);
1167                             break;
1168                         }
1169                         case 4: // lsdisc (discrete)
1170                         {
1171                             strMeth = "lsdisc";
1172                             C2F(lsdisc)(ode_f, &sizeYc, pdYData, &tleft, &tright, rwork, &rworkSize, &istate);
1173                             break;
1174                         }
1175                         case 5: // lsrgk (rk)
1176                         {
1177                             strMeth = "lsrgk";
1178                             C2F(lsrgk)(ode_f, &sizeYc, pdYData, &tleft, &tright, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1179                             break;
1180                         }
1181                         case 6: // rkf45 (rkf)
1182                         {
1183                             strMeth = "rkf45";
1184                             C2F(rkf45)(ode_f, &sizeYc, pdYData, &tleft, &tright, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1185                             break;
1186                         }
1187                         case 7: // rksimp (fix)
1188                         {
1189                             strMeth = "rksimp";
1190                             C2F(rksimp)(ode_f, &sizeYc, pdYData, &tleft, &tright, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1191                             break;
1192                         }
1193                         default: // case 0: // lsoda
1194                         {
1195                             strMeth = "lsoda";
1196                             C2F(lsoda)(ode_f, &sizeYc, pdYData, &tleft, &tright, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jt);
1197                         }
1198                     }
1199
1200                     // check error
1201                     err = checkOdeError(meth, istate);
1202                     if (err == 1) // error case
1203                     {
1204                         Scierror(999, _("%s: %s exit with state %d.\n"), "odedc", strMeth.c_str(), istate);
1205                     }
1206                 }
1207                 catch (ast::ScilabError &e)
1208                 {
1209                     char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
1210                     sciprint(_("%s: exception caught in '%s' subroutine.\n"), "odedc", strMeth.c_str());
1211                     Scierror(999, pstrMsg);
1212                     FREE(pstrMsg);
1213                     err = 1;
1214                 }
1215
1216                 if (err == 1) // error case
1217                 {
1218                     DifferentialEquation::removeDifferentialEquationFunctions();
1219                     FREE(pdYData);
1220                     FREE(YSize);
1221                     if (jroot)
1222                     {
1223                         FREE(jroot);
1224                     }
1225                     if (dStructTab)
1226                     {
1227                         FREE(dStructTab);
1228                     }
1229                     if (iStructTab)
1230                     {
1231                         FREE(iStructTab);
1232                     }
1233                     if (itol == 1 || itol == 3)
1234                     {
1235                         FREE(atol);
1236                     }
1237                     if (itol < 3)
1238                     {
1239                         FREE(rtol);
1240                     }
1241                     return types::Function::Error;
1242                 }
1243
1244                 if (err == 2) // warning case
1245                 {
1246                     if (getWarningMode())
1247                     {
1248                         sciprint(_("Integration was stoped at t = %lf.\n"), tleft);
1249                     }
1250                     break;
1251                 }
1252
1253                 double* copy = (double*)MALLOC(*YSize * sizeof(double));
1254                 C2F(dcopy)(YSize, pdYData, &one, copy, &one);
1255                 pDblYOutList.push_back(copy);
1256                 pDblTOutList.push_back(tleft);
1257
1258                 if (meth == 3 && istate == 3 && getWarningMode())
1259                 {
1260                     // istate == 3  means the integration was successful, and one or more
1261                     //              roots were found before satisfying the stop condition
1262                     //              specified by itask.  see jroot.
1263
1264                     sciprint(_("%s: Warning: At t = %lf, y is a root, jroot = "), "odedc", tleft);
1265                     for (int k = 0; k < pDblNg->get(0); k++)
1266                     {
1267                         sciprint("\t%d", jroot[k]);
1268                     }
1269                     sciprint("\n");
1270                     break;
1271                 }
1272             }
1273         }
1274         while (fabs(tleft - pDblT->get(iLastT)) > 1.0e-12 || bIntegrateContPart);
1275
1276         // Create result
1277         pDblYOut = new types::Double(pDblY0->getRows() + 1, (int)pDblYOutList.size() * pDblY0->getCols());
1278         int iSizeList = (int)pDblYOutList.size();
1279         for (int i = 0; i < iSizeList; i++)
1280         {
1281             pDblYOut->set(i * (*YSize + 1), pDblTOutList.front());
1282             for (int j = 0; j < *YSize; j++)
1283             {
1284                 pDblYOut->set(i * (*YSize + 1) + (j + 1), pDblYOutList.front()[j]);
1285             }
1286             pDblYOutList.pop_front();
1287             pDblTOutList.pop_front();
1288         }
1289     }
1290     else
1291     {
1292         // Create result
1293         pDblYOut = new types::Double(pDblY0->getRows(), pDblT->getSize() * pDblY0->getCols());
1294
1295         for (int i = 0; i < pDblT->getSize(); i++)
1296         {
1297             hf = pDblT0->get(0) + (nhpass + delta) * pDblStdel->get(0);
1298             tf = pDblT->get(i);
1299
1300             if (pDblOdeOptions && pDblOdeOptions->get(9) == 1)
1301             {
1302                 sciprint("tf - hf = %lf\n", tf - hf);
1303             }
1304
1305             if (fabs(tf - hf) < 1.0e-12)
1306             {
1307                 tright  = hf;
1308                 bStore  = true;
1309                 bUpdate = true;
1310                 if (pDblOdeOptions && pDblOdeOptions->get(9) == 1)
1311                 {
1312                     sciprint(_("integ. from tleft= %lf to hf=tf= %lf\n"), tleft, tright);
1313                 }
1314             }
1315             else if (tf < hf)
1316             {
1317                 tright  = tf;
1318                 bStore  = true;
1319                 bUpdate = false;
1320                 if (pDblOdeOptions && pDblOdeOptions->get(9) == 1)
1321                 {
1322                     sciprint(_("integ. from tleft= %lf to tf= %lf\n"), tleft, tright);
1323                 }
1324             }
1325             else if (tf > hf)
1326             {
1327                 tright  = hf;
1328                 bStore  = false;
1329                 bUpdate = true;
1330                 if (pDblOdeOptions && pDblOdeOptions->get(9) == 1)
1331                 {
1332                     sciprint(_("integ. from tleft= %lf to hf= %lf\n"), tleft, tright);
1333                 }
1334             }
1335
1336             rwork[0] = hf; // tcrit = hf
1337             try
1338             {
1339                 switch (meth)
1340                 {
1341                     case 1: // lsode (adams)
1342                     case 2: // lsode (stiff)
1343                     {
1344                         strMeth = "lsode";
1345                         int jacType = 10 * meth + jt;
1346                         C2F(lsode)(ode_f, &sizeYc, pdYData, &tleft, &tright, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jacType);
1347                         break;
1348                     }
1349                     case 3: // lsodar
1350                     {
1351                         strMeth = "lsodar";
1352                         int ng = (int)pDblNg->get(0);
1353                         C2F(lsodar)(ode_f, &sizeYc, pdYData, &tleft, &tright, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jt, ode_g, &ng, jroot);
1354                         break;
1355                     }
1356                     case 4: // lsdisc (discrete)
1357                     {
1358                         strMeth = "lsdisc";
1359                         C2F(lsdisc)(ode_f, &sizeYc, pdYData, &tleft, &tright, rwork, &rworkSize, &istate);
1360                         break;
1361                     }
1362                     case 5: // lsrgk (rk)
1363                     {
1364                         strMeth = "lsrgk";
1365                         C2F(lsrgk)(ode_f, &sizeYc, pdYData, &tleft, &tright, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1366                         break;
1367                     }
1368                     case 6: // rkf45 (rkf)
1369                     {
1370                         strMeth = "rkf45";
1371                         C2F(rkf45)(ode_f, &sizeYc, pdYData, &tleft, &tright, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1372                         break;
1373                     }
1374                     case 7: // rksimp (fix)
1375                     {
1376                         strMeth = "rksimp";
1377                         C2F(rksimp)(ode_f, &sizeYc, pdYData, &tleft, &tright, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1378                         break;
1379                     }
1380                     default: // case 0: // lsoda
1381                     {
1382                         strMeth = "lsoda";
1383                         C2F(lsoda)(ode_f, &sizeYc, pdYData, &tleft, &tright, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jt);
1384                         break;
1385                     }
1386                 }
1387
1388                 // check error
1389                 err = checkOdeError(meth, istate);
1390                 if (err == 1) // error case
1391                 {
1392                     Scierror(999, _("%s: %s exit with state %d.\n"), "odedc", strMeth.c_str(), istate);
1393                 }
1394             }
1395             catch (ast::ScilabError &e)
1396             {
1397                 char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
1398                 sciprint(_("%s: exception caught in '%s' subroutine.\n"), "odedc", strMeth.c_str());
1399                 Scierror(999, pstrMsg);
1400                 FREE(pstrMsg);
1401                 err = 1;
1402             }
1403
1404             if (err == 1) // error case
1405             {
1406                 DifferentialEquation::removeDifferentialEquationFunctions();
1407                 FREE(pdYData);
1408                 FREE(YSize);
1409                 if (jroot)
1410                 {
1411                     FREE(jroot);
1412                 }
1413                 if (dStructTab)
1414                 {
1415                     FREE(dStructTab);
1416                 }
1417                 if (iStructTab)
1418                 {
1419                     FREE(iStructTab);
1420                 }
1421                 if (itol == 1 || itol == 3)
1422                 {
1423                     FREE(atol);
1424                 }
1425                 if (itol < 3)
1426                 {
1427                     FREE(rtol);
1428                 }
1429                 return types::Function::Error;
1430             }
1431
1432             if (bUpdate)
1433             {
1434                 deFunctionsManager->setOdedcFlag();
1435
1436                 if (pDblOdeOptions && pDblOdeOptions->get(9) == 1)
1437                 {
1438                     sciprint(_("update at t = %lf\n"), tright);
1439                 }
1440
1441                 try
1442                 {
1443                     ode_f(&sizeYc, &tright, pdYData, pdYData + sizeYc);
1444                 }
1445                 catch (ast::ScilabError &e)
1446                 {
1447                     char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
1448                     sciprint(_("%s: Update failed at t = %lf\n"), "odedc", tright);
1449                     Scierror(999, pstrMsg);
1450                     DifferentialEquation::removeDifferentialEquationFunctions();
1451                     FREE(pdYData);
1452                     FREE(YSize);
1453                     if (jroot)
1454                     {
1455                         FREE(jroot);
1456                     }
1457                     if (dStructTab)
1458                     {
1459                         FREE(dStructTab);
1460                     }
1461                     if (iStructTab)
1462                     {
1463                         FREE(iStructTab);
1464                     }
1465                     if (itol == 1 || itol == 3)
1466                     {
1467                         FREE(atol);
1468                     }
1469                     if (itol < 3)
1470                     {
1471                         FREE(rtol);
1472                     }
1473                     return types::Function::Error;
1474                 }
1475
1476                 deFunctionsManager->resetOdedcFlag();
1477                 nhpass++;
1478             }
1479
1480             if (bStore)
1481             {
1482                 for (int j = 0; j < *YSize; j++)
1483                 {
1484                     pDblYOut->set(i * (*YSize) + j, pdYData[j]);
1485                 }
1486             }
1487             else
1488             {
1489                 i--;
1490             }
1491
1492             if (err == 2) // warning case
1493             {
1494                 if (getWarningMode())
1495                 {
1496                     sciprint(_("Integration was stoped at t = %lf.\n"), tleft);
1497                 }
1498
1499                 types::Double* pDblYOutTemp = pDblYOut;
1500                 pDblYOut = new types::Double(pDblYOutTemp->getRows(), i + 1);
1501                 for (int k = 0; k < pDblYOut->getSize(); k++)
1502                 {
1503                     pDblYOut->set(k, pDblYOutTemp->get(k));
1504                 }
1505                 delete pDblYOutTemp;
1506                 break;
1507             }
1508
1509             if (meth == 3 && istate == 3 && getWarningMode())
1510             {
1511                 // istate == 3  means the integration was successful, and one or more
1512                 //              roots were found before satisfying the stop condition
1513                 //              specified by itask.  see jroot.
1514
1515                 sciprint(_("%s: Warning: At t = %lf, y is a root, jroot = "), "odedc", tleft);
1516                 for (int k = 0; k < pDblNg->get(0); k++)
1517                 {
1518                     sciprint("\t%d", jroot[k]);
1519                 }
1520                 sciprint("\n");
1521
1522                 types::Double* pDblYOutTemp = pDblYOut;
1523                 pDblYOut = new types::Double(pDblYOutTemp->getRows(), i + 1);
1524                 for (int k = 0; k < pDblYOut->getSize(); k++)
1525                 {
1526                     pDblYOut->set(k, pDblYOutTemp->get(k));
1527                 }
1528                 delete pDblYOutTemp;
1529                 break;
1530             }
1531         }
1532     }
1533
1534     if (meth < 4)
1535     {
1536         if (_iRetCount > 2) //save ls0001 and eh0001 following w and iw.
1537         {
1538             int dPos    = 0;
1539             int iPos    = 0;
1540             int dSize   = 219;
1541
1542             if (dStructTab == NULL)
1543             {
1544                 dStructTab = (double*)MALLOC(dStructTabSize * sizeof(double));
1545                 iStructTab = (int*)MALLOC(iStructTabSize * sizeof(int));
1546             }
1547
1548             C2F(dcopy)(&dSize, ls0001d, &one, dStructTab, &one);
1549             memcpy(iStructTab, ls0001i, 39 * sizeof(int));
1550
1551             dPos = dSize;
1552             iPos = 39;
1553
1554             //save lsa001
1555             if (meth == 0 || meth == 3)
1556             {
1557                 dSize = 22;
1558                 C2F(dcopy)(&dSize, lsa001d, &one, dStructTab + dPos, &one);
1559                 memcpy(&iStructTab[iPos], lsa001i, 9 * sizeof(int));
1560
1561                 dPos += dSize;
1562                 iPos += 9;
1563             }
1564
1565             //save lsr001
1566             if (meth == 3)
1567             {
1568                 dSize = 5;
1569                 C2F(dcopy)(&dSize, lsr001d, &one, dStructTab + dPos, &one);
1570                 memcpy(&iStructTab[iPos], lsr001i, 9 * sizeof(int));
1571
1572                 iPos += 9;
1573             }
1574
1575             memcpy(&iStructTab[iPos], eh0001i, 2 * sizeof(int));
1576         }
1577     }
1578
1579     // *** Return result in Scilab. ***
1580
1581     out.push_back(pDblYOut); // y
1582
1583     if (meth == 3 && _iRetCount >= 2)
1584     {
1585         int sizeOfRd = 1;
1586         int k = 0;
1587
1588         for (int i = 0; i < pDblNg->get(0); i++)
1589         {
1590             if (jroot[i])
1591             {
1592                 sizeOfRd++;
1593             }
1594         }
1595
1596         types::Double* pDblRd = new types::Double(1, sizeOfRd);
1597         //rd: The first entry contains the stopping time.
1598         pDblRd->set(0, C2F(lsr001).tlast);
1599         for (int i = 0; i < pDblNg->get(0); i++)
1600         {
1601             if (jroot[i])
1602             {
1603                 k++;
1604                 pDblRd->set(k, (double)i + 1);
1605             }
1606         }
1607         out.push_back(pDblRd); // rd
1608     }
1609
1610     if ((meth < 3 && _iRetCount == 3) || (meth == 3 && _iRetCount == 4))
1611     {
1612         types::Double* pDblWOut = new types::Double(1, rwSize);
1613         C2F(dcopy)(&rworkSize, rwork, &one, pDblWOut->get(), &one);
1614         C2F(dcopy)(&dStructTabSize, dStructTab, &one, pDblWOut->get() + rworkSize, &one);
1615
1616         types::Double* pDblIwOut = new types::Double(1, iwSize);
1617         for (int i = 0; i < iworkSize; i++)
1618         {
1619             pDblIwOut->set(i, (double)iwork[i]);
1620         }
1621
1622         for (int i = 0; i < iStructTabSize; i++)
1623         {
1624             pDblIwOut->set(iworkSize + i, (double)iStructTab[i]);
1625         }
1626
1627         out.push_back(pDblWOut);  // w
1628         out.push_back(pDblIwOut); // iw
1629     }
1630
1631     // *** FREE. ***
1632     if (itol == 1 || itol == 3) // atol is scalar
1633     {
1634         FREE(atol);
1635     }
1636
1637     if (itol < 3) // rtol is scalar
1638     {
1639         FREE(rtol);
1640     }
1641
1642     FREE(pdYData);
1643     FREE(YSize);
1644     FREE(rwork);
1645     FREE(iwork);
1646
1647     if (jroot)
1648     {
1649         FREE(jroot);
1650     }
1651
1652     if (dStructTab)
1653     {
1654         FREE(dStructTab);
1655         FREE(iStructTab);
1656     }
1657
1658     DifferentialEquation::removeDifferentialEquationFunctions();
1659
1660     return types::Function::OK;
1661 }
1662