Delphi Fortran DLL передача процедуры в качестве параметра

[gimlis]

Уважаемые Гуру, очень нужна Ваша помощь.

Стоит задача переделать программу на Фортране в подпрограмму, получить DLL и далее обращаться к ней из Дельфи.

Получилась работающая версия:

      subroutine radau5solver(neqn,t,y,rtol,atol,h0) 
      !DEC$ ATTRIBUTES STDCALL, DLLEXPORT :: radau5solver 
      !DEC$ ATTRIBUTES ALIAS : "_radau5solver" :: radau5solver 
      !DEC$ ATTRIBUTES VALUE :: neqn 
      !DEC$ ATTRIBUTES REFERENCE :: t
      !DEC$ ATTRIBUTES REFERENCE :: y
      !DEC$ ATTRIBUTES REFERENCE :: rtol
      !DEC$ ATTRIBUTES REFERENCE :: atol
      !DEC$ ATTRIBUTES VALUE :: h0

      integer neqn
      double precision y(md),dy(md),t(0:md), 
     +                 h0,rtol(md),atol(md)

type TDArrayR = array of Real;
...
procedure radau5solver(
                      neqn: Integer;
                      t: TDArrayR;
                      y: TDArrayR;
                      rtol: TDArrayR;
                      atol: TDArrayR;
                      h0: Real
                      ); stdcall; external 'radau5solver.dll' name 'radau5solver';
...
radau5solver(neqn,t,y,rtol,atol,h0);

В бывшей программе на Фортране есть основная расчетная подпрограмма:

            call radau5(neqn,odef,t(i),y,t(i+1),h, 
     +                  rtol,atol,itol, 
     +                  odejac ,ijac,mljac,mujac, 
     +                  odemas ,imas,mlmas,mumas, 
     +                  solout,iout, 
     +                  work,lwork,iwork,liwork,rpar,ipar,idid)

      subroutine odef(n,x,y,f,rpar,ipar) 
      integer n,ipar(*) 
      double precision x,y(n),f(n),rpar(*) 
      integer ierr 
      ierr = 0 
      call feval(n,x,y,y,f,ierr,rpar,ipar) 
      if (ierr.ne.0) then 
         print *, 'RADAU5D: ERROR: ', 
     +            'RADAU5 can not handle FEVAL IERR' 
         stop 
      endif 
      return 
      end 

      subroutine feval(neqn,t,y,yprime,f,ierr,rpar,ipar)
      integer neqn,ierr,ipar(*)
      double precision t,y(neqn),yprime(neqn),f(neqn),rpar(*)

      f(1) = y(2)
      f(2) = ((1-y(1)**2)*y(2)-y(1))/1.0d-6
      
      return
      end

Эти функции представляют собой систему дифференциальных уравнений. В моей основной задаче будут меняться и вид уравнений(от задачи к задаче), и коэффициенты(внутри одной задачи много раз). Оставлять реализацию в Фортране не хотелось бы.
Есть ли какие-нибудь идеи, как этого добиться?

Моя попытка выглядит следующим образом. Пробуем подменить подпрограмму odef из Дельфи.

      interface
      subroutine odef(n,x,y,f,rpar,ipar) 
      integer n,ipar(*) 
      double precision x,y(n),f(n),rpar(*) 
      end subroutine odef
      end interface

c      subroutine odef(n,x,y,f,rpar,ipar) 
c      integer n,ipar(*) 
c      double precision x,y(n),f(n),rpar(*) 
c      integer ierr 
c      ierr = 0 
c      call feval(n,x,y,y,f,ierr,rpar,ipar) 
c      if (ierr.ne.0) then 
c         print *, 'RADAU5D: ERROR: ', 
c     +            'RADAU5 can not handle FEVAL IERR' 
c         stop 
c      endif 
c      return 
c      end

type TDArrayR = array of Real;
     Todef = procedure (
              n: Integer;
              x: Real;
              y: TDArrayR;
              f: TDArrayR;
              rpar: array of Real;
              ipar: array of Integer
              ); stdcall;

procedure radau5solver(
                      neqn: Integer;
                      t: TDArrayR;
                      y: TDArrayR;
                      odef: Todef;
                      rtol: TDArrayR;
                      atol: TDArrayR;
                      h0: Real
                      ); stdcall; external 'radau5solver.dll' name 'radau5solver';

procedure odef(
              n: Integer;
              x: Real;
              y: TDArrayR;
              f: TDArrayR;
              rpar: array of Real;
              ipar: array of Integer
              ); stdcall;
begin                                  //сюда возвращает, ошибка при переходе к следующему шагу
  ShowMessage('OoOoo');
//  f[0] := y[1];
//  f[1] := ((1-sqr(y[0]))*y[1]-y[0])/0.000001;
end;

radau5solver(neqn,t,y,odef,rtol,atol,h0);

При запуске программы получается следующее: запускается подпрограмма, подгружается DLL, выполняются все строчки DLL до call odef(...), далее возврат в Дельфи на begin procedure odef, а на следующем шаге выскакивает ошибка EAccessViolation "Access violation at address XXX in module 'YYY'. Read of address 00000011".
Подскажите, пожалуйста, в чем причина и как заставить работать.
Почему ошибка при переходе к следующему шагу после входа в процедуру?

[gimlis]

Так тоже не работает:

      interface
      subroutine odef(n,x,y,f,rpar,ipar) 
      !DEC$ ATTRIBUTES STDCALL, DLLEXPORT :: odef 
      !DEC$ ATTRIBUTES ALIAS : "_odef" :: odef 
      !DEC$ ATTRIBUTES VALUE :: n 
      !DEC$ ATTRIBUTES VALUE :: x
      !DEC$ ATTRIBUTES REFERENCE :: y
      !DEC$ ATTRIBUTES REFERENCE :: f
      !DEC$ ATTRIBUTES REFERENCE :: rpar
      !DEC$ ATTRIBUTES REFERENCE :: ipar       
      integer n,ipar(*) 
      double precision x,y(n),f(n),rpar(*) 
      end subroutine odef
      end interface