# Slight change in an equation brings Domain error with mcpSolve

**URL:** <https://discourse.julialang.org/t/slight-change-in-an-equation-brings-domain-error-with-mcpsolve/13881>\
**Category:** Modelling & Simulations\
**Created:** [August 22, 2018, 1:04pm UTC](https://discourse.julialang.org/t/slight-change-in-an-equation-brings-domain-error-with-mcpsolve/13881 "2018-08-22T13:04:51Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![axelle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/axelle/32/4488_2.png) [@axelle](https://discourse.julialang.org/u/axelle)\
**Post date:** [August 22, 2018, 1:04pm UTC](https://discourse.julialang.org/t/slight-change-in-an-equation-brings-domain-error-with-mcpsolve/13881/1 "2018-08-22T13:04:51Z")

</div>

Good afternoon,

I have a system of equations and keep getting “domain error” for some values of the parameters. I have been able to solve a lot of them. For instance, I had some difficulties with some integral that were not defined for values of unknow near zero that I (probably artificially) solved by imposing the unknowns to be above 0.01 instead of zero.

But the latest domain error below, I really do not understand.

When i change equation F[1] of the sytem from `F[1]= x[2] * D `, it works i.e. I do not get domain error (- the algorithm does not converge but as this is a simplified version of my true system, I am not concerned about that for now).

But when I change it to `F[1]= x[2] * D + 2 `, I get a domain error (so this is even before looking into convergence for a solution).

Would anyone know why? Understanding what is going on here would help me solve this error and the future ones. Many thanks for your help!

```julia
using QuantEcon
using NLsolve
using ForwardDiff

Re=5.0
Rl=2.0
E=4.0
D=0.1

function foc2(θlow::AbstractFloat,
                D::AbstractFloat) 
    
v(c) = 4 * c
c21(θ) = 1/ (1 - θ)   
val2(θ) = θ * D * v(c21(θ))
val22 = quadgk(val2, 0, θlow)[1]

return val22 
    
end 

 function f!(F ,x) 
            F[1]= x[2] * D
            F[2]= foc2(x[9], D)
            F[3]= D - x[3] - x[1] 
            F[4]= E - x[4] - D 
            F[5]= x[2] * D - x[5] * ((Re-Rl)/Rl) 
            F[6]= x[4] - x[6] - x[5]
            F[7]= x[8] - x[5] / x[1]  
            F[8]= x[9] - x[3]/ (D * x[2]) 
            F[9]= x[2] * D * x[7] + x[2] * D * (1 - x[7]) * (x[8]/Re) - x[3] - x[1] * x[8] 
       
end 
        
        
    initial_x = [D; D; 0.0; E/2.0; 0.01; E/2.0; 0.3; 0.1; 0.1] 
    initial_F = similar(initial_x)
        
    df = OnceDifferentiable(f!, initial_x, initial_F)
       
    results1_bk = mcpsolve(df, [0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01 , 0.0], 
                        [D, D, D, E, E, D, 0.9, 0.99, 0.9], initial_x )
 

```