diff --git a/PartitionedSolvers/CHANGELOG.md b/PartitionedSolvers/CHANGELOG.md
index 98bdf322..caa8c134 100644
--- a/PartitionedSolvers/CHANGELOG.md
+++ b/PartitionedSolvers/CHANGELOG.md
@@ -6,6 +6,13 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [0.3.1] - 2024-11-11
+
+### Added
+
+- Function `default_solver`
+- First draft of a Newton-Raphson solver in function `newton_raphson`.
+
 ## [0.3.0] - 2024-10-29
 
 ### Changed
diff --git a/PartitionedSolvers/Project.toml b/PartitionedSolvers/Project.toml
index 7a6a21dd..a420d123 100644
--- a/PartitionedSolvers/Project.toml
+++ b/PartitionedSolvers/Project.toml
@@ -1,7 +1,7 @@
 name = "PartitionedSolvers"
 uuid = "11b65f7f-80ac-401b-9ef2-3db765482d62"
 authors = ["Francesc Verdugo <f.verdugo.rojano@vu.nl>"]
-version = "0.3.0"
+version = "0.3.1"
 
 [deps]
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
diff --git a/PartitionedSolvers/src/PartitionedSolvers.jl b/PartitionedSolvers/src/PartitionedSolvers.jl
index 5ea3a80e..20c3ebe6 100644
--- a/PartitionedSolvers/src/PartitionedSolvers.jl
+++ b/PartitionedSolvers/src/PartitionedSolvers.jl
@@ -13,5 +13,6 @@ include("interfaces.jl")
 include("wrappers.jl")
 include("smoothers.jl")
 include("amg.jl")
+include("nonlinear_solvers.jl")
 
 end # module
diff --git a/PartitionedSolvers/src/interfaces.jl b/PartitionedSolvers/src/interfaces.jl
index fca4998e..27eb9dd6 100644
--- a/PartitionedSolvers/src/interfaces.jl
+++ b/PartitionedSolvers/src/interfaces.jl
@@ -26,7 +26,7 @@ abstract type AbstractSolver <: AbstractType end
 #    end
 #end
 
-function solve(solver;kwargs...)
+function solve(solver::AbstractSolver;kwargs...)
     solver, state = step(solver;kwargs...)
     while state !== :stop
         solver, state = step(solver,state)
@@ -34,14 +34,29 @@ function solve(solver;kwargs...)
     solver
 end
 
-function history(solver;kwargs...)
+function history(solver::AbstractSolver;kwargs...)
     History(identity,solver,kwargs)
 end
 
-function history(f,solver;kwargs...)
+function history(f,solver::AbstractSolver;kwargs...)
     History(f,solver,kwargs)
 end
 
+function solve(p::AbstractProblem;kwargs...)
+    s = default_solver(p)
+    solve(s;kwargs...)
+end
+
+function history(p::AbstractProblem;kwargs...)
+    s = default_solver(p)
+    history(s;kwargs...)
+end
+
+function history(f,p::AbstractProblem;kwargs...)
+    s = default_solver(p)
+    history(f,s;kwargs...)
+end
+
 struct History{A,B,C}
     f::A
     solver::B
@@ -95,6 +110,8 @@ jacobian(a::AbstractSolver) = jacobian(problem(a))
 
 abstract type AbstractLinearProblem <: AbstractProblem end
 
+default_solver(p::AbstractLinearProblem) = LinearAlgebra_lu(p)
+
 #struct LinearProblemAge <: AbstractAge
 #    solution::Int
 #    matrix::Int
@@ -229,6 +246,8 @@ end
 
 abstract type AbstractNonlinearProblem <: AbstractProblem end
 
+default_solver(p::AbstractNonlinearProblem) = newton_raphson(p)
+
 #function update(p::AbstractNonlinearProblem;kwargs...)
 #    function update_nonlinear_problem(x)
 #        update(p,x;kwargs...)
diff --git a/PartitionedSolvers/src/nonlinear_solvers.jl b/PartitionedSolvers/src/nonlinear_solvers.jl
new file mode 100644
index 00000000..48eb3e14
--- /dev/null
+++ b/PartitionedSolvers/src/nonlinear_solvers.jl
@@ -0,0 +1,78 @@
+# TODO This is a very vanilla NR solver for the moment
+function newton_raphson(p;
+        solver=default_solver,
+        iterations=1000,
+        residual_tol = convert(real(eltype(solution(p))), 1e-8),
+        solution_tol = zero(real(eltype(solution(p)))),
+        verbose = false,
+        output_prefix = "",
+    )
+    @assert uses_mutable_types(p)
+    iteration = 0
+    x = solution(p)
+    # Just to make sure that the problem is in the right state
+    p = update(p,solution=x)
+    dx = similar(x,axes(jacobian(p),2))
+    lp = linear_problem(dx,jacobian(p),residual(p))
+    ls = solver(lp)
+    solution_loss = maximum(abs,dx)
+    residual_loss = maximum(abs,residual(p))
+    workspace = (;ls,iteration,iterations,solution_loss,residual_loss,solution_tol,residual_tol,verbose,output_prefix)
+    nonlinear_solver(
+        newton_raphson_update,
+        newton_raphson_step,
+        p,
+        workspace)
+end
+
+function newton_raphson_update(ws,p)
+    (;ls,iteration,iterations,solution_loss,residual_loss,solution_tol,residual_tol,verbose,output_prefix) = ws
+    ls = update(ls;matrix=jacobian(p),rhs=residual(p))
+    (;ls,iteration,iterations,solution_loss,residual_loss,solution_tol,residual_tol,verbose,output_prefix)
+end
+
+function newton_raphson_step(ws,p,phase=:start)
+    (;ls,iteration,iterations,solution_loss,residual_loss,solution_tol,residual_tol,verbose,output_prefix) = ws
+    if phase === :start
+        iteration = 0
+        ws = (;ls,iteration,iterations,solution_loss,residual_loss,solution_tol,residual_tol,verbose,output_prefix)
+        phase = :advance
+        print_progress_header(ws)
+        print_progress(ws)
+    end
+    ls = solve(ls)
+    dx = solution(ls)
+    x = solution(p)
+    x .-= dx
+    p = update(p,solution=x)
+    iteration += 1
+    solution_loss = maximum(abs,dx)
+    residual_loss = maximum(abs,residual(p))
+    if solution_loss <= solution_tol || residual_loss <= residual_tol
+        phase = :stop
+    end
+    if iteration == iterations
+        phase = :stop
+    end
+    ws = (;ls,iteration,iterations,solution_loss,residual_loss,solution_tol,residual_tol,verbose,output_prefix)
+    ws = newton_raphson_update(ws,p)
+    print_progress(ws)
+    ws,p,phase
+end
+
+function print_progress_header(a)
+    s = a.output_prefix
+    v = a.verbose
+    c = "current"
+    t = "target"
+    v && @printf "%s%20s %20s %20s\n" s "iterations" "residual" "solution"
+    v && @printf "%s%10s %10s %10s %10s %10s %10s\n" s c t c t c t
+end
+
+function print_progress(a)
+    s = a.output_prefix
+    v = a.verbose
+    v && @printf "%s%10i %10i %10.2e %10.2e %10.2e %10.2e\n" s a.iteration a.iterations a.residual_loss a.residual_tol a.solution_loss a.solution_tol
+end
+
+
diff --git a/PartitionedSolvers/test/nonlinear_solvers_tests.jl b/PartitionedSolvers/test/nonlinear_solvers_tests.jl
new file mode 100644
index 00000000..801607c7
--- /dev/null
+++ b/PartitionedSolvers/test/nonlinear_solvers_tests.jl
@@ -0,0 +1,49 @@
+module NonlinearSolversTests
+
+import PartitionedSolvers as PS
+import PartitionedArrays as PA
+using Test
+using LinearAlgebra
+
+function mock_nonlinear_problem(x0)
+    r0 = zeros(1)
+    j0 = PA.sparse_matrix([1],[1],[0.0],1,1)
+    workspace = nothing
+    PS.nonlinear_problem(x0,r0,j0,workspace) do p
+        x = PS.solution(p)
+        r = PS.residual(p)
+        j = PS.jacobian(p)
+        if r !== nothing
+            r .= 2 .* x.^2 .- 4
+            p = PS.update(p,residual = r)
+        end
+        if j !== nothing
+            j .= 4 .* x
+            p = PS.update(p,jacobian = j)
+        end
+        p
+    end |> PS.update
+end
+
+x = [1.0]
+p = mock_nonlinear_problem(x)
+s = PS.newton_raphson(p;verbose=true)
+s = PS.solve(s)
+
+x = [1.0]
+p = mock_nonlinear_problem(x)
+s = PS.solve(p)
+
+x = [1.0]
+p = mock_nonlinear_problem(x)
+for s in PS.history(p)
+    @show s.workspace.residual_loss
+end
+
+x = [1.0]
+p = mock_nonlinear_problem(x)
+for residual_loss in PS.history(s->s.workspace.residual_loss,p)
+    @show residual_loss
+end
+
+end # module
diff --git a/PartitionedSolvers/test/runtests.jl b/PartitionedSolvers/test/runtests.jl
index a4ffbe83..3da585dd 100644
--- a/PartitionedSolvers/test/runtests.jl
+++ b/PartitionedSolvers/test/runtests.jl
@@ -7,6 +7,7 @@ using Test
 @testset "PartitionedSolvers" begin
     @testset "interfaces" begin include("interfaces_tests.jl") end
     @testset "wrappers" begin include("wrappers_tests.jl") end
+    @testset "nonlinear_solvers" begin include("nonlinear_solvers_tests.jl") end
     @testset "smoothers" begin include("smoothers_tests.jl") end
     @testset "amg" begin include("amg_tests.jl") end
 end
diff --git a/PartitionedSolvers/test/smoothers_tests.jl b/PartitionedSolvers/test/smoothers_tests.jl
index 4a34c673..a9179013 100644
--- a/PartitionedSolvers/test/smoothers_tests.jl
+++ b/PartitionedSolvers/test/smoothers_tests.jl
@@ -15,6 +15,14 @@ A = psparse(args...) |> fetch
 x = pones(partition(axes(A,2)))
 b = A*x
 
+tol = 1.e-8
+y = similar(x)
+y .= 0
+p = PS.linear_problem(y,A,b)
+s = PS.solve(p)
+y = PS.solution(s)
+@test norm(y-x)/norm(x) < tol
+
 tol = 1.e-8
 y = similar(x)
 y .= 0