# Using SciMLOperators in LinearSolve

**URL:** https://discourse.julialang.org/t/using-scimloperators-in-linearsolve/104596
**Category:** Numerics
**Tags:** question
**Created:** [October 4, 2023, 11:52pm UTC](https://discourse.julialang.org/t/using-scimloperators-in-linearsolve/104596 "2023-10-04T23:52:32Z")
**Posts on this page:** 1
**Page:** 1

<div class="post-metadata">

### Author: ![j35](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j35/32/202256_2.png) [@j35](https://discourse.julialang.org/u/j35)
#### Post date: [October 4, 2023, 11:52pm UTC](https://discourse.julialang.org/t/using-scimloperators-in-linearsolve/104596/1 "2023-10-04T23:52:32Z")

</div>

I’m trying to use FunctionOperators from SciMLOperators, to solve a linear problem using matrix free method. But I’m getting wrong result. What am I doing wrong here? Thanks in advance.

Here is the MWE.

```julia
using SparseArrays
using LinearAlgebra
using LinearSolve
using SciMLOperators

# du = A * (B * u)
f(u, p, t) = p[1] * (p[2] * u)
F = FunctionOperator(f, zeros(n), zeros(n))

n = 100
u = rand(n)
u0 = zeros(n)
A, B = rand(n, n), rand(n, n)
ref = B \ (A \ u)
p = (A, B)

relative_residue_ref = norm(u - (A * (B * ref))) / norm(u) # 6.910495879970626e-14
relative_residue_function = norm(u - f(ref, p, Float64[])) / norm(u) # 6.910495879970626e-14

prob = LinearProblem(F, u, p=p)
linsolve = init(prob)

sol1 = solve(linsolve)
relative_residue_function_operator = norm(u - (A * (B * sol1.u))) / norm(u) # 384.0793068887865

```