If2 used in PDE solution in simultaneous dynamic mode

I am currently working on a model consists of a system of partial differential equations (PDE). The differentials that vary in space have been discretized, while differentials in time have been solved numerically by Orthogonal Collocation on Finite Elements method.

I have a problem. When using condition if2, the initial condition is zero, regardless of what was informed.

It is necessary to initialize the model with IMODE=4 and then with IMODE = 5, both with SOLVER = 3.

This is a simplified example of the model:

import numpy as np

from gekko import GEKKO
m = GEKKO(remote=False)

zf = 0.17
tf = 1

n = 3
nt = 2
m.time = np.linspace(0,tf,nt)

dz = zf/n

#########################

tg_ini = [580.8821, 326.10596036, 318.07166567]
tg = [m.Var(tg_ini[i]) for i in range(n)]

cpv2_ini = [2.0038420499, 1.872134158, 1.8696804845]
cpv2 = [m.Param(cpv2_ini[i]) for i in range(n)]

tp = [m.Var(310.15) for i in range(n)]

#########################

def func_cpv11(t):   
    cp11 = ( 31.32234 + ((-20.23531)*(t/1000)) + (57.86644*((t/1000)**2)) + ((-36.50624)*((t/1000)**3)) + ((-0.007374)/((t/1000)**2)) )/31.998
    return cp11     

def func_cpv21(t):
    cp21 = ( 28.98641 + (1.853978*(t/1000)) + ((-9.647459)*((t/1000)**2)) + ((16.63537)*((t/1000)**3)) + (0.000117/((t/1000)**2)) )/28.0134
    return cp21

def func_cpv22(t): 
    cp22 = ( 19.50583 + (19.88705*(t/1000)) + ((-8.598535)*((t/1000)**2)) + (1.369784*((t/1000)**3)) + (0.527601/((t/1000)**2)) )/28.0134
    return cp22

#########################

cpv1 = [m.Intermediate(func_cpv11(tg[i])) for i in range(n)]   

cpv21 = [m.Intermediate(func_cpv21(tg[i])) for i in range(n)]
cpv22 = [m.Intermediate(func_cpv22(tg[i])) for i in range(n)]

for i in range(n):
    cpv2[i] = m.if2(tg[i]-500,cpv21[i],cpv22[i])

cpv = [[cpv1[i], cpv2[i]] for i in range(n)]

y = np.array([0.5, 0.5])

fg_cpv = [m.Intermediate(y @ cpv[i]) for i in range(n)]

#########################

m.Equation(tg[0] == 580.8821)
m.Equation([tg[i] == tg[i-1] + dz*fg_cpv[i-1] for i in range(1,n)])

m.Equation([tp[i].dt() == fg_cpv[i] for i in range(n)])

#########################

m.options.imode = 5
m.options.SOLVER = 3

m.solve()

print('cpv2 = ') 
for i in range(n):
    print(cpv2[i].value)

This is the result obtained:

----------------------------------------------------------------
 APMonitor, Version 0.9.2
 APMonitor Optimization Suite
 ----------------------------------------------------------------
 
 
 --------- APM Model Size ------------
 Each time step contains
   Objects      :  3
   Constants    :  0
   Variables    :  21
   Intermediates:  12
   Connections  :  6
   Equations    :  27
   Residuals    :  15
 
 Number of state variables:    27
 Number of total equations: -  21
 Number of slack variables: -  0
 ---------------------------------------
 Degrees of freedom       :    6
 
 **********************************************
 Dynamic Estimation with Interior Point Solver
 **********************************************
  
  
 Info: Exact Hessian

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

This is Ipopt version 3.10.2, running with linear solver mumps.

Number of nonzeros in equality constraint Jacobian...:       52
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:       12

Total number of variables............................:       27
                     variables with only lower bounds:        6
                variables with lower and upper bounds:       12
                     variables with only upper bounds:        0
Total number of equality constraints.................:       21
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0 6.0000000e+003 2.55e+002 9.90e+001   0.0 0.00e+000    -  0.00e+000 0.00e+000   0
   1 9.2649845e+003 2.55e+002 9.90e+001  11.2 2.64e+011    -  8.34e-012 3.76e-012f  1
   2 3.6289712e+004 2.55e+002 5.95e+004  12.5 3.20e+012    -  9.33e-010 1.41e-012f  1
   3 3.1244934e+006 2.55e+002 3.55e+011  11.8 1.21e+010    -  1.00e+000 4.26e-008f  1
   4 1.2233739e+008 6.90e-002 7.62e+011  11.1 1.99e+004    -  5.67e-001 1.00e+000f  1
   5 1.2233421e+008 2.16e-012 6.44e+009   5.4 5.33e-001    -  9.90e-001 1.00e+000f  1
   6 1.2198510e+008 2.63e-012 6.44e+007   3.4 5.82e+001    -  9.90e-001 1.00e+000f  1
   7 9.4816241e+007 1.14e-012 6.46e+005   1.4 4.53e+003    -  9.90e-001 1.00e+000f  1
   8 3.0450882e+006 1.21e-012 8.97e+003  -0.3 1.53e+004    -  9.86e-001 1.00e+000f  1
   9 1.6038743e+005 3.18e-011 4.47e+002   2.8 5.28e+002    -  9.96e-001 9.03e-001f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10 3.5917803e+004 2.72e-011 3.42e+002   2.6 3.23e+001    -  1.00e+000 7.43e-001f  1
  11 4.4032685e+003 3.28e-011 2.67e-001   1.5 1.15e+000    -  9.99e-001 1.00e+000f  1
  12 4.9020198e+001 7.06e-013 7.78e-003  -0.5 1.87e-001    -  9.94e-001 9.95e-001f  1
  13 4.8932248e-001 1.14e-013 8.96e-005  -6.3 1.94e-003    -  9.90e-001 9.90e-001f  1
  14 -1.7447063e-003 3.93e-014 9.21e-002  -8.3 1.94e-005    -  9.99e-001 9.96e-001f  1
  15 -2.2492197e-003 1.14e-013 3.63e+003  -5.9 7.80e-008    -  1.00e+000 5.42e-001f  1
  16 -2.4686751e-003 1.14e-013 2.28e+002  -6.7 4.09e-008    -  1.00e+000 9.58e-001f  1
  17 -2.4879150e-003 1.14e-013 2.92e-001  -8.3 2.71e-009    -  1.00e+000 9.99e-001f  1
  18 -2.4882162e-003 1.14e-013 9.09e-013 -10.1 2.70e-011    -  1.00e+000 1.00e+000h  1

Number of Iterations....: 18

                                   (scaled)                 (unscaled)
Objective...............: -2.4882162366554839e-004  -2.4882162366554835e-003
Dual infeasibility......:  9.0949470177292824e-013   9.0949470177292824e-012
Constraint violation....:  3.4861931869050499e-014   1.1368683772161603e-013
Complementarity.........:  7.5538837469831659e-011   7.5538837469831659e-010
Overall NLP error.......:  9.1075642737176940e-012   7.5538837469831659e-010


Number of objective function evaluations             = 19
Number of objective gradient evaluations             = 19
Number of equality constraint evaluations            = 19
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 19
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 18
Total CPU secs in IPOPT (w/o function evaluations)   =      0.053
Total CPU secs in NLP function evaluations           =      0.020

EXIT: Optimal Solution Found.

 The solution was found.

 The final value of the objective function is  -0.0024882162366554835
 
 ---------------------------------------------------
 Solver         :  IPOPT (v3.12)
 Solution time  :  0.10469999999999999 sec
 Objective      :  0.
 Successful solution
 ---------------------------------------------------
 

cpv2 = 
[0.0, 1.0705093271]
[0.0, 1.0705216696]
[0.0, 1.0705340141]

It is a simple model and I cannot understand the cause of the problem.