diff --git a/thys/Matrices_for_ODEs/MTX_Examples.thy b/thys/Matrices_for_ODEs/MTX_Examples.thy
new file mode 100644
--- /dev/null
+++ b/thys/Matrices_for_ODEs/MTX_Examples.thy
@@ -0,0 +1,237 @@
+(*  Title:       Verification Examples
+    Author:      Jonathan Julián Huerta y Munive, 2020
+    Maintainer:  Jonathan Julián Huerta y Munive <jjhuertaymunive1@sheffield.ac.uk>
+*)
+
+section \<open> Verification examples \<close>
+
+
+theory MTX_Examples
+  imports MTX_Flows "Hybrid_Systems_VCs.HS_VC_Spartan"
+
+begin
+
+
+subsection \<open> Examples \<close>
+
+abbreviation hoareT :: "('a \<Rightarrow> bool) \<Rightarrow> ('a \<Rightarrow> 'a set) \<Rightarrow> ('a \<Rightarrow> bool) \<Rightarrow> bool" 
+  ("PRE_ HP _ POST _" [85,85]85) where "PRE P HP X POST Q \<equiv> (P \<le> |X]Q)"
+
+
+subsubsection \<open> Verification by uniqueness. \<close>
+
+abbreviation mtx_circ :: "2 sq_mtx" ("A")
+  where "A \<equiv> mtx  
+   ([0,  1] # 
+    [-1, 0] # [])"
+
+abbreviation mtx_circ_flow :: "real \<Rightarrow> real^2 \<Rightarrow> real^2" ("\<phi>")
+  where "\<phi> t s \<equiv> (\<chi> i. if i = 1 then s$1 * cos t + s$2 * sin t else - s$1 * sin t + s$2 * cos t)"
+
+lemma mtx_circ_flow_eq: "exp (t *\<^sub>R A) *\<^sub>V s = \<phi> t s"
+  apply(rule local_flow.eq_solution[OF local_flow_sq_mtx_linear, symmetric])
+    apply(rule ivp_solsI, simp_all add: sq_mtx_vec_mult_eq vec_eq_iff)
+  unfolding UNIV_2 using exhaust_2
+  by (force intro!: poly_derivatives simp: matrix_vector_mult_def)+
+
+lemma mtx_circ: 
+  "PRE(\<lambda>s. r\<^sup>2 = (s $ 1)\<^sup>2 + (s $ 2)\<^sup>2) 
+  HP x\<acute>=(*\<^sub>V) A & G 
+  POST (\<lambda>s. r\<^sup>2 = (s $ 1)\<^sup>2 + (s $ 2)\<^sup>2)"
+  apply(subst local_flow.fbox_g_ode[OF local_flow_sq_mtx_linear])
+  unfolding mtx_circ_flow_eq by auto
+
+no_notation mtx_circ ("A")
+        and mtx_circ_flow ("\<phi>")
+
+
+subsubsection \<open> Flow of diagonalisable matrix. \<close>
+
+abbreviation mtx_hOsc :: "real \<Rightarrow> real \<Rightarrow> 2 sq_mtx" ("A")
+  where "A a b \<equiv> mtx  
+   ([0, 1] # 
+    [a, b] # [])"
+
+abbreviation mtx_chB_hOsc :: "real \<Rightarrow> real \<Rightarrow> 2 sq_mtx" ("P")
+  where "P a b \<equiv> mtx
+   ([a, b] # 
+    [1, 1] # [])"
+
+lemma inv_mtx_chB_hOsc: 
+  "a \<noteq> b \<Longrightarrow> (P a b)\<^sup>-\<^sup>1 = (1/(a - b)) *\<^sub>R mtx 
+   ([ 1, -b] # 
+    [-1,  a] # [])"
+  apply(rule sq_mtx_inv_unique, unfold scaleR_mtx2 times_mtx2)
+  by (simp add: diff_divide_distrib[symmetric] one_mtx2)+
+
+lemma invertible_mtx_chB_hOsc: "a \<noteq> b \<Longrightarrow> mtx_invertible (P a b)"
+  apply(rule mtx_invertibleI[of _ "(P a b)\<^sup>-\<^sup>1"])
+   apply(unfold inv_mtx_chB_hOsc scaleR_mtx2 times_mtx2 one_mtx2)
+  by (subst sq_mtx_eq_iff, simp add: vector_def frac_diff_eq1)+
+
+lemma mtx_hOsc_diagonalizable:
+  fixes a b :: real
+  defines "\<iota>\<^sub>1 \<equiv> (b - sqrt (b^2+4*a))/2" and "\<iota>\<^sub>2 \<equiv> (b + sqrt (b^2+4*a))/2"
+  assumes "b\<^sup>2 + a * 4 > 0" and "a \<noteq> 0"
+  shows "A a b = P (-\<iota>\<^sub>2/a) (-\<iota>\<^sub>1/a) * (\<d>\<i>\<a>\<g> i. if i = 1 then \<iota>\<^sub>1 else \<iota>\<^sub>2) * (P (-\<iota>\<^sub>2/a) (-\<iota>\<^sub>1/a))\<^sup>-\<^sup>1"
+  unfolding assms apply(subst inv_mtx_chB_hOsc)
+  using assms(3,4) apply(simp_all add: diag2_eq[symmetric])
+  unfolding sq_mtx_times_eq sq_mtx_scaleR_eq UNIV_2 apply(subst sq_mtx_eq_iff)
+  using exhaust_2 assms by (auto simp: field_simps, auto simp: field_power_simps)
+
+lemma mtx_hOsc_solution_eq:
+  fixes a b :: real
+  defines "\<iota>\<^sub>1 \<equiv> (b - sqrt (b\<^sup>2+4*a))/2" and "\<iota>\<^sub>2 \<equiv> (b + sqrt (b\<^sup>2+4*a))/2"
+  defines "\<Phi> t \<equiv> mtx (
+   [\<iota>\<^sub>2*exp(t*\<iota>\<^sub>1) - \<iota>\<^sub>1*exp(t*\<iota>\<^sub>2),     exp(t*\<iota>\<^sub>2)-exp(t*\<iota>\<^sub>1)]#
+   [a*exp(t*\<iota>\<^sub>2) - a*exp(t*\<iota>\<^sub>1), \<iota>\<^sub>2*exp(t*\<iota>\<^sub>2)-\<iota>\<^sub>1*exp(t*\<iota>\<^sub>1)]#[])"
+  assumes "b\<^sup>2 + a * 4 > 0" and "a \<noteq> 0"
+  shows "P (-\<iota>\<^sub>2/a) (-\<iota>\<^sub>1/a) * (\<d>\<i>\<a>\<g> i. exp (t * (if i=1 then \<iota>\<^sub>1 else \<iota>\<^sub>2))) * (P (-\<iota>\<^sub>2/a) (-\<iota>\<^sub>1/a))\<^sup>-\<^sup>1 
+  = (1/sqrt (b\<^sup>2 + a * 4)) *\<^sub>R (\<Phi> t)"
+  unfolding assms apply(subst inv_mtx_chB_hOsc)
+  using assms apply(simp_all add: mtx_times_scaleR_commute, subst sq_mtx_eq_iff)
+  unfolding UNIV_2 sq_mtx_times_eq sq_mtx_scaleR_eq sq_mtx_uminus_eq apply(simp_all add: axis_def)
+  by (auto simp: field_simps, auto simp: field_power_simps)+
+ 
+lemma local_flow_mtx_hOsc:
+  fixes a b
+  defines "\<iota>\<^sub>1 \<equiv> (b - sqrt (b^2+4*a))/2" and "\<iota>\<^sub>2 \<equiv> (b + sqrt (b^2+4*a))/2"
+  defines "\<Phi> t \<equiv> mtx (
+   [\<iota>\<^sub>2*exp(t*\<iota>\<^sub>1) - \<iota>\<^sub>1*exp(t*\<iota>\<^sub>2),     exp(t*\<iota>\<^sub>2)-exp(t*\<iota>\<^sub>1)]#
+   [a*exp(t*\<iota>\<^sub>2) - a*exp(t*\<iota>\<^sub>1), \<iota>\<^sub>2*exp(t*\<iota>\<^sub>2)-\<iota>\<^sub>1*exp(t*\<iota>\<^sub>1)]#[])"
+  assumes "b\<^sup>2 + a * 4 > 0" and "a \<noteq> 0"
+  shows "local_flow ((*\<^sub>V) (A a b)) UNIV UNIV (\<lambda>t. (*\<^sub>V) ((1/sqrt (b\<^sup>2 + a * 4)) *\<^sub>R \<Phi> t))"
+  unfolding assms using local_flow_sq_mtx_linear[of "A a b"] assms
+  apply(subst (asm) exp_scaleR_diagonal2[OF invertible_mtx_chB_hOsc mtx_hOsc_diagonalizable])
+     apply(simp, simp, simp)
+  by (subst (asm) mtx_hOsc_solution_eq) simp_all
+
+lemma overdamped_door_arith:
+  assumes "b\<^sup>2 + a * 4 > 0" and "a < 0" and "b \<le> 0" and "t \<ge> 0" and "s1 > 0"
+  shows "0 \<le> ((b + sqrt (b\<^sup>2 + 4 * a)) * exp (t * (b - sqrt (b\<^sup>2 + 4 * a)) / 2) / 2 - 
+(b - sqrt (b\<^sup>2 + 4 * a)) * exp (t * (b + sqrt (b\<^sup>2 + 4 * a)) / 2) / 2) * s1 / sqrt (b\<^sup>2 + a * 4)"
+proof(subst diff_divide_distrib[symmetric], simp)
+  have f0: "s1 / (2 * sqrt (b\<^sup>2 + a * 4)) > 0"  (is "s1/?c3 > 0")
+    using assms(1,5) by simp
+  have f1: "(b - sqrt (b\<^sup>2 + 4 * a)) < (b + sqrt (b\<^sup>2 + 4 * a))" (is "?c2 < ?c1") 
+    and f2: "(b + sqrt (b\<^sup>2 + 4 * a)) < 0"
+    using sqrt_ge_absD[of b "b\<^sup>2 + 4 * a"] assms by (force, linarith)
+  hence f3: "exp (t * ?c2 / 2) \<le> exp (t * ?c1 / 2)" (is "exp ?t1 \<le> exp ?t2")
+    unfolding exp_le_cancel_iff 
+    using assms(4) by (case_tac "t=0", simp_all)
+  hence "?c2 * exp ?t2 \<le> ?c2 * exp ?t1"
+    using f1 f2 real_mult_le_cancel_iff2[of "-?c2" "exp ?t1" "exp ?t2"] by linarith 
+  also have "... < ?c1 * exp ?t1"
+    using f1 by auto
+  also have"... \<le> ?c1 * exp ?t1"
+    using f1 f2 by auto
+  ultimately show "0 \<le> (?c1 * exp ?t1 - ?c2 * exp ?t2) * s1 / ?c3"
+    using f0 f1 assms(5) by auto
+qed
+
+lemma overdamped_door:
+  assumes "b\<^sup>2 + a * 4 > 0" and "a < 0" and "b \<le> 0" and "0 \<le> t"
+  shows "PRE (\<lambda>s. s$1 = 0)
+  HP (LOOP 
+      (\<lambda>s. {s. s$1 > 0 \<and> s$2 = 0});
+      (x\<acute>=(*\<^sub>V) (A a b) & G on {0..t} UNIV @ 0) 
+     INV (\<lambda>s. 0 \<le> s$1))
+  POST (\<lambda>s. 0 \<le> s $ 1)"
+  apply(rule fbox_loopI, simp_all add: le_fun_def)
+  apply(subst local_flow.fbox_g_ode_ivl[OF local_flow_mtx_hOsc[OF assms(1)]])
+  using assms apply(simp_all add: le_fun_def fbox_def)
+  unfolding sq_mtx_scaleR_eq UNIV_2 sq_mtx_vec_mult_eq
+  by (clarsimp simp: overdamped_door_arith)
+
+
+no_notation mtx_hOsc ("A")
+        and mtx_chB_hOsc ("P")
+
+
+subsubsection \<open> Flow of non-diagonalisable matrix. \<close>
+
+abbreviation mtx_cnst_acc :: "3 sq_mtx" ("K")
+  where "K \<equiv> mtx (
+  [0,1,0] #
+  [0,0,1] # 
+  [0,0,0] # [])"
+
+lemma pow2_scaleR_mtx_cnst_acc: "(t *\<^sub>R K)\<^sup>2 = mtx (
+  [0,0,t\<^sup>2] #
+  [0,0,0] # 
+  [0,0,0] # [])"
+  unfolding power2_eq_square apply(subst sq_mtx_eq_iff)
+  unfolding sq_mtx_times_eq UNIV_3 by auto
+
+lemma powN_scaleR_mtx_cnst_acc: "n > 2 \<Longrightarrow> (t *\<^sub>R K)^n = 0"
+  apply(induct n, simp, case_tac "n \<le> 2")
+   apply(subgoal_tac "n = 2", erule ssubst)
+  unfolding power_Suc2 pow2_scaleR_mtx_cnst_acc sq_mtx_times_eq UNIV_3
+  by (auto simp: sq_mtx_eq_iff)
+
+lemma exp_mtx_cnst_acc: "exp (t *\<^sub>R K) = ((t *\<^sub>R K)\<^sup>2/\<^sub>R 2) + (t *\<^sub>R K) + 1"
+  unfolding exp_def apply(subst suminf_eq_sum[of 2])
+  using powN_scaleR_mtx_cnst_acc by (simp_all add: numeral_2_eq_2)
+ 
+lemma exp_mtx_cnst_acc_simps:
+  "exp (t *\<^sub>R K) $$ 1 $ 1 = 1" "exp (t *\<^sub>R K) $$ 1 $ 2 = t" "exp (t *\<^sub>R K) $$ 1 $ 3 = t^2/2"
+  "exp (t *\<^sub>R K) $$ 2 $ 1 = 0" "exp (t *\<^sub>R K) $$ 2 $ 2 = 1" "exp (t *\<^sub>R K) $$ 2 $ 3 = t"
+  "exp (t *\<^sub>R K) $$ 3 $ 1 = 0" "exp (t *\<^sub>R K) $$ 3 $ 2 = 0" "exp (t *\<^sub>R K) $$ 3 $ 3 = 1"
+  unfolding exp_mtx_cnst_acc one_mtx3 pow2_scaleR_mtx_cnst_acc by simp_all
+
+lemma exp_mtx_cnst_acc_vec_mult_eq: "exp (t *\<^sub>R K) *\<^sub>V s = 
+  vector [s$3 * t^2/2 + s$2 * t + s$1, s$3 * t + s$2, s$3]"
+  apply(simp add: sq_mtx_vec_mult_eq vector_def)
+  unfolding UNIV_3 apply (simp add: exp_mtx_cnst_acc_simps fun_eq_iff)
+  using exhaust_3 exp_mtx_cnst_acc_simps(7,8,9) by fastforce
+
+lemma local_flow_mtx_cnst_acc:
+  "local_flow ((*\<^sub>V) K) UNIV UNIV (\<lambda>t s. ((t *\<^sub>R K)\<^sup>2/\<^sub>R 2 + (t *\<^sub>R K) + 1) *\<^sub>V s)"
+  using local_flow_sq_mtx_linear[of K] unfolding exp_mtx_cnst_acc .  
+
+lemma docking_station_arith:
+  assumes "(d::real) > x" and "v > 0"
+  shows "(v = v\<^sup>2 * t / (2 * d - 2 * x)) \<longleftrightarrow> (v * t - v\<^sup>2 * t\<^sup>2 / (4 * d - 4 * x) + x = d)"
+proof
+  assume "v = v\<^sup>2 * t / (2 * d - 2 * x)"
+  hence "v * t = 2 * (d - x)"
+    using assms by (simp add: eq_divide_eq power2_eq_square) 
+  hence "v * t - v\<^sup>2 * t\<^sup>2 / (4 * d - 4 * x) + x = 2 * (d - x) - 4 * (d - x)\<^sup>2 / (4 * (d - x)) + x"
+    apply(subst power_mult_distrib[symmetric])
+    by (erule ssubst, subst power_mult_distrib, simp)
+  also have "... = d"
+    apply(simp only: mult_divide_mult_cancel_left_if)
+    using assms by (auto simp: power2_eq_square)
+  finally show "v * t - v\<^sup>2 * t\<^sup>2 / (4 * d - 4 * x) + x = d" .
+next
+  assume "v * t - v\<^sup>2 * t\<^sup>2 / (4 * d - 4 * x) + x = d"
+  hence "0 = v\<^sup>2 * t\<^sup>2 / (4 * (d - x)) + (d - x) - v * t"
+    by auto
+  hence "0 = (4 * (d - x)) * (v\<^sup>2 * t\<^sup>2 / (4 * (d - x)) + (d - x) - v * t)"
+    by auto
+  also have "... = v\<^sup>2 * t\<^sup>2 + 4 * (d - x)\<^sup>2  - (4 * (d - x)) * (v * t)"
+    using assms apply(simp add: distrib_left right_diff_distrib)
+    apply(subst right_diff_distrib[symmetric])+
+    by (simp add: power2_eq_square)
+  also have "... = (v * t - 2 * (d - x))\<^sup>2"
+    by (simp only: power2_diff, auto simp: field_simps power2_diff)
+  finally have "0 = (v * t - 2 * (d - x))\<^sup>2" .
+  hence "v * t = 2 * (d - x)"
+    by auto
+  thus "v = v\<^sup>2 * t / (2 * d - 2 * x)"
+    apply(subst power2_eq_square, subst mult.assoc)
+    apply(erule ssubst, subst right_diff_distrib[symmetric])
+    using assms by auto
+qed
+
+lemma docking_station:
+  assumes "d > x\<^sub>0" and "v\<^sub>0 > 0"
+  shows "PRE (\<lambda>s. s$1 = x\<^sub>0 \<and> s$2 = v\<^sub>0)
+  HP ((3 ::= (\<lambda>s. -(v\<^sub>0^2/(2*(d-x\<^sub>0))))); x\<acute>=(*\<^sub>V) K & G)
+  POST (\<lambda>s. s$2 = 0 \<longleftrightarrow> s$1 = d)"
+  apply(clarsimp simp: le_fun_def local_flow.fbox_g_ode[OF local_flow_sq_mtx_linear[of K]])
+  unfolding exp_mtx_cnst_acc_vec_mult_eq using assms by (simp add: docking_station_arith)
+
+no_notation mtx_cnst_acc ("K")
+
+end
\ No newline at end of file
diff --git a/thys/Matrices_for_ODEs/MTX_Flows.thy b/thys/Matrices_for_ODEs/MTX_Flows.thy
new file mode 100644
--- /dev/null
+++ b/thys/Matrices_for_ODEs/MTX_Flows.thy
@@ -0,0 +1,291 @@
+(*  Title:       Affine systems of ODEs
+    Author:      Jonathan Julián Huerta y Munive, 2020
+    Maintainer:  Jonathan Julián Huerta y Munive <jjhuertaymunive1@sheffield.ac.uk>
+*)
+
+section \<open> Affine systems of ODEs \<close>
+
+text \<open>Affine systems of ordinary differential equations (ODEs) are those whose vector 
+fields are linear operators. Broadly speaking, if there are functions $A$ and $B$ such that the 
+system of ODEs $X'\, t = f\, (X\, t)$ turns into $X'\, t = (A\, t)\cdot (X\, t)+(B\, t)$, then it
+is affine. The end goal of this section is to prove that every affine system of ODEs has a unique 
+solution, and to obtain a characterization of said solution. \<close>
+
+theory MTX_Flows
+  imports SQ_MTX "Hybrid_Systems_VCs.HS_ODEs"
+
+begin
+
+
+subsection \<open> Existence and uniqueness for affine systems \<close>
+
+definition matrix_continuous_on :: "real set \<Rightarrow> (real \<Rightarrow> ('a::real_normed_algebra_1)^'n^'m) \<Rightarrow> bool" 
+  where "matrix_continuous_on T A = (\<forall>t \<in> T. \<forall>\<epsilon> > 0. \<exists> \<delta> > 0. \<forall>\<tau>\<in>T. \<bar>\<tau> - t\<bar> < \<delta> \<longrightarrow> \<parallel>A \<tau> - A t\<parallel>\<^sub>o\<^sub>p \<le> \<epsilon>)"
+
+lemma continuous_on_matrix_vector_multl:
+  assumes "matrix_continuous_on T A"
+  shows "continuous_on T (\<lambda>t. A t *v s)"
+proof(rule continuous_onI, simp add: dist_norm)
+  fix e t::real assume "0 < e" and "t \<in> T"
+  let ?\<epsilon> = "e/(\<parallel>(if s = 0 then 1 else s)\<parallel>)"
+  have "?\<epsilon> > 0"
+    using \<open>0 < e\<close> by simp 
+  then obtain \<delta> where dHyp: "\<delta> > 0 \<and> (\<forall>\<tau>\<in>T. \<bar>\<tau> - t\<bar> < \<delta> \<longrightarrow> \<parallel>A \<tau> - A t\<parallel>\<^sub>o\<^sub>p \<le> ?\<epsilon>)"
+    using assms \<open>t \<in> T\<close> unfolding dist_norm matrix_continuous_on_def by fastforce
+  {fix \<tau> assume "\<tau> \<in> T" and "\<bar>\<tau> - t\<bar> < \<delta>"
+    have obs: "?\<epsilon> * (\<parallel>s\<parallel>) = (if s = 0 then 0 else e)"
+      by auto
+    have "\<parallel>A \<tau> *v s - A t *v s\<parallel> = \<parallel>(A \<tau> - A t) *v s\<parallel>"
+      by (simp add: matrix_vector_mult_diff_rdistrib)      
+    also have "... \<le> (\<parallel>A \<tau> - A t\<parallel>\<^sub>o\<^sub>p) * (\<parallel>s\<parallel>)"
+      using norm_matrix_le_mult_op_norm by blast
+    also have "... \<le> ?\<epsilon> * (\<parallel>s\<parallel>)"
+      using dHyp \<open>\<tau> \<in> T\<close> \<open>\<bar>\<tau> - t\<bar> < \<delta>\<close> mult_right_mono norm_ge_zero by blast 
+    finally have "\<parallel>A \<tau> *v s - A t *v s\<parallel> \<le> e"
+      by (subst (asm) obs) (metis (mono_tags, hide_lams) \<open>0 < e\<close> less_eq_real_def order_trans)}
+  thus "\<exists>d>0. \<forall>\<tau>\<in>T. \<bar>\<tau> - t\<bar> < d \<longrightarrow> \<parallel>A \<tau> *v s - A t *v s\<parallel> \<le> e"
+    using dHyp by blast
+qed
+
+lemma lipschitz_cond_affine:
+  fixes A :: "real \<Rightarrow> 'a::real_normed_algebra_1^'n^'m" and T::"real set"
+  defines "L \<equiv> Sup {\<parallel>A t\<parallel>\<^sub>o\<^sub>p |t. t \<in> T}"
+  assumes "t \<in> T" and "bdd_above {\<parallel>A t\<parallel>\<^sub>o\<^sub>p |t. t \<in> T}"
+  shows "\<parallel>A t *v x - A t *v y\<parallel> \<le> L * (\<parallel>x - y\<parallel>)"
+proof-
+  have obs: "\<parallel>A t\<parallel>\<^sub>o\<^sub>p \<le> Sup {\<parallel>A t\<parallel>\<^sub>o\<^sub>p |t. t \<in> T}"
+    apply(rule cSup_upper)
+    using continuous_on_subset assms by (auto simp: dist_norm)
+  have "\<parallel>A t *v x - A t *v y\<parallel> = \<parallel>A t *v (x - y)\<parallel>"
+    by (simp add: matrix_vector_mult_diff_distrib)
+  also have "... \<le> (\<parallel>A t\<parallel>\<^sub>o\<^sub>p) * (\<parallel>x - y\<parallel>)"
+    using norm_matrix_le_mult_op_norm by blast
+  also have "... \<le> Sup {\<parallel>A t\<parallel>\<^sub>o\<^sub>p |t. t \<in> T} * (\<parallel>x - y\<parallel>)"
+    using obs mult_right_mono norm_ge_zero by blast 
+  finally show "\<parallel>A t *v x - A t *v y\<parallel> \<le> L * (\<parallel>x - y\<parallel>)"
+    unfolding assms .
+qed
+
+lemma local_lipschitz_affine:
+  fixes A :: "real \<Rightarrow> 'a::real_normed_algebra_1^'n^'m"
+  assumes "open T" and "open S" 
+    and Ahyp: "\<And>\<tau> \<epsilon>. \<epsilon> > 0 \<Longrightarrow> \<tau> \<in> T \<Longrightarrow> cball \<tau> \<epsilon> \<subseteq> T \<Longrightarrow> bdd_above {\<parallel>A t\<parallel>\<^sub>o\<^sub>p |t. t \<in> cball \<tau> \<epsilon>}"
+  shows "local_lipschitz T S (\<lambda>t s. A t *v s + B t)"
+proof(unfold local_lipschitz_def lipschitz_on_def, clarsimp)
+  fix s t assume "s \<in> S" and "t \<in> T"
+  then obtain e1 e2 where "cball t e1 \<subseteq> T" and "cball s e2 \<subseteq> S" and "min e1 e2 > 0"
+    using open_cballE[OF _ \<open>open T\<close>] open_cballE[OF _ \<open>open S\<close>] by force
+  hence obs: "cball t (min e1 e2) \<subseteq> T"
+    by auto
+  let ?L = "Sup {\<parallel>A \<tau>\<parallel>\<^sub>o\<^sub>p |\<tau>. \<tau> \<in> cball t (min e1 e2)}"
+  have "\<parallel>A t\<parallel>\<^sub>o\<^sub>p \<in> {\<parallel>A \<tau>\<parallel>\<^sub>o\<^sub>p |\<tau>. \<tau> \<in> cball t (min e1 e2)}"
+    using \<open>min e1 e2 > 0\<close> by auto
+  moreover have bdd: "bdd_above {\<parallel>A \<tau>\<parallel>\<^sub>o\<^sub>p |\<tau>. \<tau> \<in> cball t (min e1 e2)}"
+    by (rule Ahyp, simp only: \<open>min e1 e2 > 0\<close>, simp_all add: \<open>t \<in> T\<close> obs)
+  moreover have "Sup {\<parallel>A \<tau>\<parallel>\<^sub>o\<^sub>p |\<tau>. \<tau> \<in> cball t (min e1 e2)} \<ge> 0"
+    apply(rule order.trans[OF op_norm_ge_0[of "A t"]])
+    by (rule cSup_upper[OF calculation])
+  moreover have "\<forall>x\<in>cball s (min e1 e2) \<inter> S. \<forall>y\<in>cball s (min e1 e2) \<inter> S. 
+    \<forall>\<tau>\<in>cball t (min e1 e2) \<inter> T. dist (A \<tau> *v x) (A \<tau> *v y) \<le> ?L * dist x y"
+    apply(clarify, simp only: dist_norm, rule lipschitz_cond_affine)
+    using \<open>min e1 e2 > 0\<close> bdd by auto
+  ultimately show "\<exists>e>0. \<exists>L. \<forall>t\<in>cball t e \<inter> T. 0 \<le> L \<and> 
+    (\<forall>x\<in>cball s e \<inter> S. \<forall>y\<in>cball s e \<inter> S. dist (A t *v x) (A t *v y) \<le> L * dist x y)"
+    using \<open>min e1 e2 > 0\<close> by blast
+qed
+
+lemma picard_lindeloef_affine:
+  fixes A :: "real \<Rightarrow> 'a::{banach,real_normed_algebra_1,heine_borel}^'n^'n"
+  assumes Ahyp: "matrix_continuous_on T A"
+      and "\<And>\<tau> \<epsilon>. \<tau> \<in> T \<Longrightarrow> \<epsilon> > 0 \<Longrightarrow> bdd_above {\<parallel>A t\<parallel>\<^sub>o\<^sub>p |t. dist \<tau> t \<le> \<epsilon>}"
+      and Bhyp: "continuous_on T B" and "open S" 
+      and "t\<^sub>0 \<in> T" and Thyp: "open T" "is_interval T" 
+    shows "picard_lindeloef (\<lambda> t s. A t *v s + B t) T S t\<^sub>0"
+  apply (unfold_locales, simp_all add: assms, clarsimp)
+   apply (rule continuous_on_add[OF continuous_on_matrix_vector_multl[OF Ahyp] Bhyp])
+  by (rule local_lipschitz_affine) (simp_all add: assms)
+
+lemma picard_lindeloef_autonomous_affine: 
+  fixes A :: "'a::{banach,real_normed_field,heine_borel}^'n^'n"
+  shows "picard_lindeloef (\<lambda> t s. A *v s + B) UNIV UNIV t\<^sub>0"
+  using picard_lindeloef_affine[of _ "\<lambda>t. A" "\<lambda>t. B"] 
+  unfolding matrix_continuous_on_def by (simp only: diff_self op_norm0, auto)
+
+lemma picard_lindeloef_autonomous_linear:
+  fixes A :: "'a::{banach,real_normed_field,heine_borel}^'n^'n"
+  shows "picard_lindeloef (\<lambda> t. (*v) A) UNIV UNIV t\<^sub>0"
+  using picard_lindeloef_autonomous_affine[of A 0] by force
+
+lemmas unique_sol_autonomous_affine = picard_lindeloef.unique_solution[OF 
+    picard_lindeloef_autonomous_affine _ _ funcset_UNIV UNIV_I _ _ funcset_UNIV UNIV_I]
+
+lemmas unique_sol_autonomous_linear = picard_lindeloef.unique_solution[OF 
+    picard_lindeloef_autonomous_linear _ _ funcset_UNIV UNIV_I _ _ funcset_UNIV UNIV_I]
+
+
+subsection \<open> Flow for affine systems \<close>
+
+
+subsubsection \<open> Derivative rules for square matrices \<close>
+
+lemma has_derivative_exp_scaleRl[derivative_intros]:
+  fixes f::"real \<Rightarrow> real" (* by Fabian Immler and Johannes Hölzl *)
+  assumes "D f \<mapsto> f' at t within T"
+  shows "D (\<lambda>t. exp (f t *\<^sub>R A)) \<mapsto> (\<lambda>h. f' h *\<^sub>R (exp (f t *\<^sub>R A) * A)) at t within T"
+proof -
+  have "bounded_linear f'" 
+    using assms by auto
+  then obtain m where obs: "f' = (\<lambda>h. h * m)"
+    using real_bounded_linear by blast
+  thus ?thesis
+    using vector_diff_chain_within[OF _ exp_scaleR_has_vector_derivative_right] 
+      assms obs by (auto simp: has_vector_derivative_def comp_def)
+qed
+
+lemma has_vderiv_on_exp_scaleRl:
+  assumes "D f = f' on T"
+  shows "D (\<lambda>x. exp (f x *\<^sub>R A)) = (\<lambda>x. f' x *\<^sub>R exp (f x *\<^sub>R A) * A) on T"
+  using assms unfolding has_vderiv_on_def has_vector_derivative_def apply clarsimp
+  by (rule has_derivative_exp_scaleRl, auto simp: fun_eq_iff)
+
+lemma vderiv_on_exp_scaleRlI[poly_derivatives]:
+  assumes "D f = f' on T" and "g' = (\<lambda>x. f' x *\<^sub>R exp (f x *\<^sub>R A) * A)"
+  shows "D (\<lambda>x. exp (f x *\<^sub>R A)) = g' on T"
+  using has_vderiv_on_exp_scaleRl assms by simp
+
+lemma has_derivative_mtx_ith[derivative_intros]:
+  fixes t::real and T :: "real set"
+  defines "t\<^sub>0 \<equiv> netlimit (at t within T)"
+  assumes "D A \<mapsto> (\<lambda>h. h *\<^sub>R A' t) at t within T"
+  shows "D (\<lambda>t. A t $$ i) \<mapsto> (\<lambda>h. h *\<^sub>R A' t $$ i) at t within T"
+  using assms unfolding has_derivative_def apply safe
+   apply(force simp: bounded_linear_def bounded_linear_axioms_def)
+  apply(rule_tac F="\<lambda>\<tau>. (A \<tau> - A t\<^sub>0 - (\<tau> - t\<^sub>0) *\<^sub>R A' t) /\<^sub>R (\<parallel>\<tau> - t\<^sub>0\<parallel>)" in tendsto_zero_norm_bound)
+  by (clarsimp, rule mult_left_mono, metis (no_types, lifting) norm_column_le_norm 
+      sq_mtx_minus_ith sq_mtx_scaleR_ith) simp_all
+
+lemmas has_derivative_mtx_vec_mult[derivative_intros] = 
+  bounded_bilinear.FDERIV[OF bounded_bilinear_sq_mtx_vec_mult]
+
+lemma vderiv_mtx_vec_mult_intro[poly_derivatives]:
+  assumes "D u = u' on T" and "D A = A' on T"
+      and "g = (\<lambda>t. A t *\<^sub>V u' t + A' t *\<^sub>V u t )"
+    shows "D (\<lambda>t. A t *\<^sub>V u t) = g on T"
+  using assms unfolding has_vderiv_on_def has_vector_derivative_def apply clarsimp
+  apply(erule_tac x=x in ballE, simp_all)+
+  apply(rule derivative_eq_intros(142))
+  by (auto simp: fun_eq_iff mtx_vec_scaleR_commute pth_6 scaleR_mtx_vec_assoc)
+
+lemmas has_vderiv_on_ivl_integral = ivl_integral_has_vderiv_on[OF vderiv_on_continuous_on]
+
+declare has_vderiv_on_ivl_integral [poly_derivatives]
+
+lemma has_derivative_mtx_vec_multl[derivative_intros]:
+  assumes "\<And> i j. D (\<lambda>t. (A t) $$ i $ j) \<mapsto> (\<lambda>\<tau>. \<tau> *\<^sub>R (A' t) $$ i $ j) (at t within T)"
+  shows "D (\<lambda>t. A t *\<^sub>V x) \<mapsto> (\<lambda>\<tau>. \<tau> *\<^sub>R (A' t) *\<^sub>V x) at t within T"
+  unfolding sq_mtx_vec_mult_sum_cols
+  apply(rule_tac f'1="\<lambda>i \<tau>. \<tau> *\<^sub>R  (x $ i *\<^sub>R \<c>\<o>\<l> i (A' t))" in derivative_eq_intros(10))
+   apply(simp_all add: scaleR_right.sum)
+  apply(rule_tac g'1="\<lambda>\<tau>. \<tau> *\<^sub>R \<c>\<o>\<l> i (A' t)" in derivative_eq_intros(4), simp_all add: mult.commute)
+  using assms unfolding sq_mtx_col_def column_def apply(transfer, simp)
+  apply(rule has_derivative_vec_lambda)
+  by (simp add: scaleR_vec_def)
+
+lemma continuous_on_mtx_vec_multr: "continuous_on S ((*\<^sub>V) A)"
+  by transfer (simp add: matrix_vector_mult_linear_continuous_on)
+
+\<comment> \<open>Automatically generated derivative rules from this subsubsection \<close>
+
+thm derivative_eq_intros(140,141,142,143)
+
+
+subsubsection \<open> Existence and uniqueness with square matrices \<close>
+
+text \<open>Finally, we can use the @{term exp} operation to characterize the general solutions for affine 
+systems of ODEs. We show that they satisfy the @{term "local_flow"} locale.\<close>
+
+lemma continuous_on_sq_mtx_vec_multl:
+  fixes A :: "real \<Rightarrow> ('n::finite) sq_mtx"
+  assumes "continuous_on T A"
+  shows "continuous_on T (\<lambda>t. A t *\<^sub>V s)"
+proof-
+  have "matrix_continuous_on T (\<lambda>t. to_vec (A t))"
+    using assms by (force simp: continuous_on_iff dist_norm norm_sq_mtx_def matrix_continuous_on_def)
+  hence "continuous_on T (\<lambda>t. to_vec (A t) *v s)"
+    by (rule continuous_on_matrix_vector_multl)
+  thus ?thesis
+    by transfer
+qed
+
+lemmas continuous_on_affine = continuous_on_add[OF continuous_on_sq_mtx_vec_multl]
+
+lemma local_lipschitz_sq_mtx_affine:
+  fixes A :: "real \<Rightarrow> ('n::finite) sq_mtx"
+  assumes "continuous_on T A" "open T" "open S"
+  shows "local_lipschitz T S (\<lambda>t s. A t *\<^sub>V s + B t)"
+proof-
+  have obs: "\<And>\<tau> \<epsilon>. 0 < \<epsilon> \<Longrightarrow>  \<tau> \<in> T \<Longrightarrow> cball \<tau> \<epsilon> \<subseteq> T \<Longrightarrow> bdd_above {\<parallel>A t\<parallel> |t. t \<in> cball \<tau> \<epsilon>}"
+    by (rule bdd_above_norm_cont_comp, rule continuous_on_subset[OF assms(1)], simp_all)
+  hence "\<And>\<tau> \<epsilon>. 0 < \<epsilon> \<Longrightarrow> \<tau> \<in> T \<Longrightarrow> cball \<tau> \<epsilon> \<subseteq> T \<Longrightarrow> bdd_above {\<parallel>to_vec (A t)\<parallel>\<^sub>o\<^sub>p |t. t \<in> cball \<tau> \<epsilon>}"
+    by (simp add: norm_sq_mtx_def)
+  hence "local_lipschitz T S (\<lambda>t s. to_vec (A t) *v s + B t)"
+    using local_lipschitz_affine[OF assms(2,3), of "\<lambda>t. to_vec (A t)"] by force
+  thus ?thesis
+    by transfer 
+qed
+
+lemma picard_lindeloef_sq_mtx_affine:
+  assumes "continuous_on T A" and "continuous_on T B" 
+    and "t\<^sub>0 \<in> T" "is_interval T" "open T" and "open S"
+  shows "picard_lindeloef (\<lambda>t s. A t *\<^sub>V s + B t) T S t\<^sub>0"
+  apply(unfold_locales, simp_all add: assms, clarsimp)
+  using continuous_on_affine assms apply blast
+  by (rule local_lipschitz_sq_mtx_affine, simp_all add: assms)
+
+lemmas sq_mtx_unique_sol_autonomous_affine = picard_lindeloef.unique_solution[OF 
+    picard_lindeloef_sq_mtx_affine[OF 
+      continuous_on_const 
+      continuous_on_const 
+      UNIV_I is_interval_univ 
+      open_UNIV open_UNIV] 
+    _ _ funcset_UNIV UNIV_I _ _ funcset_UNIV UNIV_I]
+
+lemma has_vderiv_on_sq_mtx_linear:
+  "D (\<lambda>t. exp ((t - t\<^sub>0) *\<^sub>R A) *\<^sub>V s) = (\<lambda>t. A *\<^sub>V (exp ((t - t\<^sub>0) *\<^sub>R A) *\<^sub>V s)) on {t\<^sub>0--t}"
+  by (rule poly_derivatives)+ (auto simp: exp_times_scaleR_commute sq_mtx_times_vec_assoc)
+
+lemma has_vderiv_on_sq_mtx_affine:
+  fixes t\<^sub>0::real and A :: "('a::finite) sq_mtx"
+  defines "lSol c t \<equiv> exp ((c * (t - t\<^sub>0)) *\<^sub>R A)"
+  shows "D (\<lambda>t. lSol 1 t *\<^sub>V s + lSol 1 t *\<^sub>V (\<integral>\<^sub>t\<^sub>0\<^sup>t (lSol (-1) \<tau> *\<^sub>V B) \<partial>\<tau>)) = 
+  (\<lambda>t. A *\<^sub>V (lSol 1 t *\<^sub>V s + lSol 1 t *\<^sub>V (\<integral>\<^sub>t\<^sub>0\<^sup>t (lSol (-1) \<tau> *\<^sub>V B) \<partial>\<tau>)) + B) on {t\<^sub>0--t}"
+  unfolding assms apply(simp only: mult.left_neutral mult_minus1)
+  apply(rule poly_derivatives, (force)?, (force)?, (force)?, (force)?)+
+  by (simp add: mtx_vec_mult_add_rdistl sq_mtx_times_vec_assoc[symmetric] 
+      exp_minus_inverse exp_times_scaleR_commute mult_exp_exp  scale_left_distrib[symmetric])
+
+lemma autonomous_linear_sol_is_exp:
+  assumes "D X = (\<lambda>t. A *\<^sub>V X t) on {t\<^sub>0--t}" and "X t\<^sub>0 = s" 
+  shows "X t = exp ((t - t\<^sub>0) *\<^sub>R A) *\<^sub>V s"
+  apply(rule sq_mtx_unique_sol_autonomous_affine[of X A 0, OF _ \<open>X t\<^sub>0 = s\<close>])
+  using assms has_vderiv_on_sq_mtx_linear by force+
+
+lemma autonomous_affine_sol_is_exp_plus_int:
+  assumes "D X = (\<lambda>t. A *\<^sub>V X t + B) on {t\<^sub>0--t}" and "X t\<^sub>0 = s" 
+  shows "X t = exp ((t - t\<^sub>0) *\<^sub>R A) *\<^sub>V s + exp ((t - t\<^sub>0) *\<^sub>R A) *\<^sub>V (\<integral>\<^sub>t\<^sub>0\<^sup>t(exp (- (\<tau> - t\<^sub>0) *\<^sub>R A) *\<^sub>V B)\<partial>\<tau>)"
+  apply(rule sq_mtx_unique_sol_autonomous_affine[OF assms])
+  using has_vderiv_on_sq_mtx_affine by force+
+
+lemma local_flow_sq_mtx_linear: "local_flow ((*\<^sub>V) A) UNIV UNIV (\<lambda>t s. exp (t *\<^sub>R A) *\<^sub>V s)"
+  unfolding local_flow_def local_flow_axioms_def apply safe
+  using picard_lindeloef_sq_mtx_affine[of _ "\<lambda>t. A" "\<lambda>t. 0"] apply force
+  using has_vderiv_on_sq_mtx_linear[of 0] by auto
+
+lemma local_flow_sq_mtx_affine: "local_flow (\<lambda>s. A *\<^sub>V s + B) UNIV UNIV 
+  (\<lambda>t s. exp (t *\<^sub>R A) *\<^sub>V s + exp (t *\<^sub>R A) *\<^sub>V (\<integral>\<^sub>0\<^sup>t(exp (- \<tau> *\<^sub>R A) *\<^sub>V B)\<partial>\<tau>))"
+  unfolding local_flow_def local_flow_axioms_def apply safe
+  using picard_lindeloef_sq_mtx_affine[of _ "\<lambda>t. A" "\<lambda>t. B"] apply force
+  using has_vderiv_on_sq_mtx_affine[of 0 A] by auto
+
+
+end
\ No newline at end of file
diff --git a/thys/Matrices_for_ODEs/MTX_Norms.thy b/thys/Matrices_for_ODEs/MTX_Norms.thy
new file mode 100644
--- /dev/null
+++ b/thys/Matrices_for_ODEs/MTX_Norms.thy
@@ -0,0 +1,285 @@
+(*  Title:       Matrix norms
+    Author:      Jonathan Julián Huerta y Munive, 2020
+    Maintainer:  Jonathan Julián Huerta y Munive <jjhuertaymunive1@sheffield.ac.uk>
+*)
+
+section \<open> Matrix norms \<close>
+
+text \<open> Here, we explore some properties about the operator and the maximum norms for matrices. \<close>
+
+theory MTX_Norms
+  imports MTX_Preliminaries
+
+begin
+
+
+subsection\<open> Matrix operator norm \<close>
+
+abbreviation op_norm :: "('a::real_normed_algebra_1)^'n^'m \<Rightarrow> real" ("(1\<parallel>_\<parallel>\<^sub>o\<^sub>p)" [65] 61)
+  where "\<parallel>A\<parallel>\<^sub>o\<^sub>p \<equiv> onorm (\<lambda>x. A *v x)"
+
+lemma norm_matrix_bound:
+  fixes A :: "('a::real_normed_algebra_1)^'n^'m"
+  shows "\<parallel>x\<parallel> = 1 \<Longrightarrow> \<parallel>A *v x\<parallel> \<le> \<parallel>(\<chi> i j. \<parallel>A $ i $ j\<parallel>) *v 1\<parallel>"
+proof-
+  fix x :: "('a, 'n) vec" assume "\<parallel>x\<parallel> = 1"
+  hence xi_le1:"\<And>i. \<parallel>x $ i\<parallel> \<le> 1" 
+    by (metis Finite_Cartesian_Product.norm_nth_le) 
+  {fix j::'m 
+    have "\<parallel>(\<Sum>i\<in>UNIV. A $ j $ i * x $ i)\<parallel> \<le> (\<Sum>i\<in>UNIV. \<parallel>A $ j $ i * x $ i\<parallel>)"
+      using norm_sum by blast
+    also have "... \<le> (\<Sum>i\<in>UNIV. (\<parallel>A $ j $ i\<parallel>) * (\<parallel>x $ i\<parallel>))"
+      by (simp add: norm_mult_ineq sum_mono)
+    also have "... \<le> (\<Sum>i\<in>UNIV. (\<parallel>A $ j $ i\<parallel>) * 1)"
+      using xi_le1 by (simp add: sum_mono mult_left_le)
+    finally have "\<parallel>(\<Sum>i\<in>UNIV. A $ j $ i * x $ i)\<parallel> \<le> (\<Sum>i\<in>UNIV. (\<parallel>A $ j $ i\<parallel>) * 1)" by simp}
+  hence "\<And>j. \<parallel>(A *v x) $ j\<parallel> \<le> ((\<chi> i1 i2. \<parallel>A $ i1 $ i2\<parallel>) *v 1) $ j"
+    unfolding matrix_vector_mult_def by simp
+  hence "(\<Sum>j\<in>UNIV. (\<parallel>(A *v x) $ j\<parallel>)\<^sup>2) \<le> (\<Sum>j\<in>UNIV. (\<parallel>((\<chi> i1 i2. \<parallel>A $ i1 $ i2\<parallel>) *v 1) $ j\<parallel>)\<^sup>2)"
+    by (metis (mono_tags, lifting) norm_ge_zero power2_abs power_mono real_norm_def sum_mono) 
+  thus "\<parallel>A *v x\<parallel> \<le> \<parallel>(\<chi> i j. \<parallel>A $ i $ j\<parallel>) *v 1\<parallel>"
+    unfolding norm_vec_def L2_set_def by simp
+qed
+
+lemma onorm_set_proptys:
+  fixes A :: "('a::real_normed_algebra_1)^'n^'m"
+  shows "bounded (range (\<lambda>x. (\<parallel>A *v x\<parallel>) / (\<parallel>x\<parallel>)))"
+    and "bdd_above (range (\<lambda>x. (\<parallel>A *v x\<parallel>) / (\<parallel>x\<parallel>)))"
+    and "(range (\<lambda>x. (\<parallel>A *v x\<parallel>) / (\<parallel>x\<parallel>))) \<noteq> {}"
+  unfolding bounded_def bdd_above_def image_def dist_real_def 
+    apply(rule_tac x=0 in exI)
+  by (rule_tac x="\<parallel>(\<chi> i j. \<parallel>A $ i $ j\<parallel>) *v 1\<parallel>" in exI, clarsimp,
+      subst mult_norm_matrix_sgn_eq[symmetric], clarsimp,
+      rule_tac x="sgn _" in norm_matrix_bound, simp add: norm_sgn)+ force
+
+lemma op_norm_set_proptys:
+  fixes A :: "('a::real_normed_algebra_1)^'n^'m"
+  shows "bounded {\<parallel>A *v x\<parallel> | x. \<parallel>x\<parallel> = 1}"
+    and "bdd_above {\<parallel>A *v x\<parallel> | x. \<parallel>x\<parallel> = 1}"
+    and "{\<parallel>A *v x\<parallel> | x. \<parallel>x\<parallel> = 1} \<noteq> {}"
+  unfolding bounded_def bdd_above_def apply safe
+    apply(rule_tac x=0 in exI, rule_tac x="\<parallel>(\<chi> i j. \<parallel>A $ i $ j\<parallel>) *v 1\<parallel>" in exI)
+    apply(force simp: norm_matrix_bound dist_real_def)
+   apply(rule_tac x="\<parallel>(\<chi> i j. \<parallel>A $ i $ j\<parallel>) *v 1\<parallel>" in exI, force simp: norm_matrix_bound)
+  using ex_norm_eq_1 by blast
+
+lemma op_norm_def: "\<parallel>A\<parallel>\<^sub>o\<^sub>p = Sup {\<parallel>A *v x\<parallel> | x. \<parallel>x\<parallel> = 1}"
+  apply(rule antisym[OF onorm_le cSup_least[OF op_norm_set_proptys(3)]])
+   apply(case_tac "x = 0", simp)
+   apply(subst mult_norm_matrix_sgn_eq[symmetric], simp)
+   apply(rule cSup_upper[OF _ op_norm_set_proptys(2)])
+   apply(force simp: norm_sgn)
+  unfolding onorm_def 
+  apply(rule cSup_upper[OF _ onorm_set_proptys(2)])
+  by (simp add: image_def, clarsimp) (metis div_by_1)
+
+lemma norm_matrix_le_op_norm: "\<parallel>x\<parallel> = 1 \<Longrightarrow> \<parallel>A *v x\<parallel> \<le> \<parallel>A\<parallel>\<^sub>o\<^sub>p"
+  apply(unfold onorm_def, rule cSup_upper[OF _ onorm_set_proptys(2)])
+  unfolding image_def by (clarsimp, rule_tac x=x in exI) simp
+
+lemma op_norm_ge_0: "0 \<le> \<parallel>A\<parallel>\<^sub>o\<^sub>p"
+  using ex_norm_eq_1 norm_ge_zero norm_matrix_le_op_norm basic_trans_rules(23) by blast
+
+lemma norm_sgn_le_op_norm: "\<parallel>A *v sgn x\<parallel> \<le> \<parallel>A\<parallel>\<^sub>o\<^sub>p"
+  by (cases "x=0", simp_all add: norm_sgn norm_matrix_le_op_norm op_norm_ge_0)
+
+lemma norm_matrix_le_mult_op_norm: "\<parallel>A *v x\<parallel> \<le> (\<parallel>A\<parallel>\<^sub>o\<^sub>p) * (\<parallel>x\<parallel>)"
+proof-
+  have "\<parallel>A *v x\<parallel> = (\<parallel>A *v sgn x\<parallel>) * (\<parallel>x\<parallel>)"
+    by(simp add: mult_norm_matrix_sgn_eq)
+  also have "... \<le> (\<parallel>A\<parallel>\<^sub>o\<^sub>p) * (\<parallel>x\<parallel>)"
+    using norm_sgn_le_op_norm[of A] by (simp add: mult_mono')
+  finally show ?thesis by simp
+qed
+
+lemma blin_matrix_vector_mult: "bounded_linear ((*v) A)" for A :: "('a::real_normed_algebra_1)^'n^'m"
+  by (unfold_locales) (auto intro: norm_matrix_le_mult_op_norm simp: 
+      mult.commute matrix_vector_right_distrib vector_scaleR_commute)
+
+lemma op_norm_eq_0: "(\<parallel>A\<parallel>\<^sub>o\<^sub>p = 0) = (A = 0)" for A :: "('a::real_normed_field)^'n^'m"
+  unfolding onorm_eq_0[OF blin_matrix_vector_mult] using matrix_axis_0[of 1 A] by fastforce
+
+lemma op_norm0: "\<parallel>(0::('a::real_normed_field)^'n^'m)\<parallel>\<^sub>o\<^sub>p = 0"
+  using op_norm_eq_0[of 0] by simp
+                  
+lemma op_norm_triangle: "\<parallel>A + B\<parallel>\<^sub>o\<^sub>p \<le> (\<parallel>A\<parallel>\<^sub>o\<^sub>p) + (\<parallel>B\<parallel>\<^sub>o\<^sub>p)" 
+  using onorm_triangle[OF blin_matrix_vector_mult[of A] blin_matrix_vector_mult[of B]]
+    matrix_vector_mult_add_rdistrib[symmetric, of A _ B] by simp
+
+lemma op_norm_scaleR: "\<parallel>c *\<^sub>R A\<parallel>\<^sub>o\<^sub>p = \<bar>c\<bar> * (\<parallel>A\<parallel>\<^sub>o\<^sub>p)"
+  unfolding onorm_scaleR[OF blin_matrix_vector_mult, symmetric] scaleR_vector_assoc ..
+
+lemma op_norm_matrix_matrix_mult_le: "\<parallel>A ** B\<parallel>\<^sub>o\<^sub>p \<le> (\<parallel>A\<parallel>\<^sub>o\<^sub>p) * (\<parallel>B\<parallel>\<^sub>o\<^sub>p)" 
+proof(rule onorm_le)
+  have "0 \<le> (\<parallel>A\<parallel>\<^sub>o\<^sub>p)" 
+    by(rule onorm_pos_le[OF blin_matrix_vector_mult])
+  fix x have "\<parallel>A ** B *v x\<parallel> = \<parallel>A *v (B *v x)\<parallel>"
+    by (simp add: matrix_vector_mul_assoc)
+  also have "... \<le> (\<parallel>A\<parallel>\<^sub>o\<^sub>p) * (\<parallel>B *v x\<parallel>)"
+    by (simp add: norm_matrix_le_mult_op_norm[of _ "B *v x"])
+  also have "... \<le> (\<parallel>A\<parallel>\<^sub>o\<^sub>p) * ((\<parallel>B\<parallel>\<^sub>o\<^sub>p) * (\<parallel>x\<parallel>))"
+    using norm_matrix_le_mult_op_norm[of B x] \<open>0 \<le> (\<parallel>A\<parallel>\<^sub>o\<^sub>p)\<close> mult_left_mono by blast
+  finally show "\<parallel>A ** B *v x\<parallel> \<le> (\<parallel>A\<parallel>\<^sub>o\<^sub>p) * (\<parallel>B\<parallel>\<^sub>o\<^sub>p) * (\<parallel>x\<parallel>)" 
+    by simp
+qed
+
+lemma norm_matrix_vec_mult_le_transpose:
+  "\<parallel>x\<parallel> = 1 \<Longrightarrow> (\<parallel>A *v x\<parallel>) \<le> sqrt (\<parallel>transpose A ** A\<parallel>\<^sub>o\<^sub>p) * (\<parallel>x\<parallel>)" for A :: "real^'n^'n"
+proof-
+  assume "\<parallel>x\<parallel> = 1"
+  have "(\<parallel>A *v x\<parallel>)\<^sup>2 = (A *v x) \<bullet> (A *v x)"
+    using dot_square_norm[of "(A *v x)"] by simp
+  also have "... = x \<bullet> (transpose A *v (A *v x))"
+    using vec_mult_inner by blast
+  also have "... \<le> (\<parallel>x\<parallel>) * (\<parallel>transpose A *v (A *v x)\<parallel>)"
+    using norm_cauchy_schwarz by blast
+  also have "... \<le> (\<parallel>transpose A ** A\<parallel>\<^sub>o\<^sub>p) * (\<parallel>x\<parallel>)^2"
+    apply(subst matrix_vector_mul_assoc) 
+    using norm_matrix_le_mult_op_norm[of "transpose A ** A" x]
+    by (simp add: \<open>\<parallel>x\<parallel> = 1\<close>) 
+  finally have "((\<parallel>A *v x\<parallel>))^2 \<le> (\<parallel>transpose A ** A\<parallel>\<^sub>o\<^sub>p) * (\<parallel>x\<parallel>)^2"
+    by linarith
+  thus "(\<parallel>A *v x\<parallel>) \<le> sqrt ((\<parallel>transpose A ** A\<parallel>\<^sub>o\<^sub>p)) * (\<parallel>x\<parallel>)"
+    by (simp add: \<open>\<parallel>x\<parallel> = 1\<close> real_le_rsqrt)
+qed
+
+lemma op_norm_le_sum_column: "\<parallel>A\<parallel>\<^sub>o\<^sub>p \<le> (\<Sum>i\<in>UNIV. \<parallel>column i A\<parallel>)" for A :: "real^'n^'m"  
+proof(unfold op_norm_def, rule cSup_least[OF op_norm_set_proptys(3)], clarsimp)
+  fix x :: "real^'n" assume x_def:"\<parallel>x\<parallel> = 1" 
+  hence x_hyp:"\<And>i. \<parallel>x $ i\<parallel> \<le> 1"
+    by (simp add: norm_bound_component_le_cart)
+  have "(\<parallel>A *v x\<parallel>) = \<parallel>(\<Sum>i\<in>UNIV. x $ i *s column i A)\<parallel>"
+    by(subst matrix_mult_sum[of A], simp)
+  also have "... \<le> (\<Sum>i\<in>UNIV. \<parallel>x $ i *s column i A\<parallel>)"
+    by (simp add: sum_norm_le)
+  also have "... = (\<Sum>i\<in>UNIV. (\<parallel>x $ i\<parallel>) * (\<parallel>column i A\<parallel>))"
+    by (simp add: mult_norm_matrix_sgn_eq)
+  also have "... \<le> (\<Sum>i\<in>UNIV. \<parallel>column i A\<parallel>)"
+    using x_hyp by (simp add: mult_left_le_one_le sum_mono) 
+  finally show "\<parallel>A *v x\<parallel> \<le> (\<Sum>i\<in>UNIV. \<parallel>column i A\<parallel>)" .
+qed
+
+lemma op_norm_le_transpose: "\<parallel>A\<parallel>\<^sub>o\<^sub>p \<le> \<parallel>transpose A\<parallel>\<^sub>o\<^sub>p" for A :: "real^'n^'n"  
+proof-
+  have obs:"\<forall>x. \<parallel>x\<parallel> = 1 \<longrightarrow> (\<parallel>A *v x\<parallel>) \<le> sqrt ((\<parallel>transpose A ** A\<parallel>\<^sub>o\<^sub>p)) * (\<parallel>x\<parallel>)"
+    using norm_matrix_vec_mult_le_transpose by blast
+  have "(\<parallel>A\<parallel>\<^sub>o\<^sub>p) \<le> sqrt ((\<parallel>transpose A ** A\<parallel>\<^sub>o\<^sub>p))"
+    using obs apply(unfold op_norm_def)
+    by (rule cSup_least[OF op_norm_set_proptys(3)]) clarsimp
+  hence "((\<parallel>A\<parallel>\<^sub>o\<^sub>p))\<^sup>2 \<le> (\<parallel>transpose A ** A\<parallel>\<^sub>o\<^sub>p)"
+    using power_mono[of "(\<parallel>A\<parallel>\<^sub>o\<^sub>p)" _ 2] op_norm_ge_0
+    by (metis not_le real_less_lsqrt)
+  also have "... \<le> (\<parallel>transpose A\<parallel>\<^sub>o\<^sub>p) * (\<parallel>A\<parallel>\<^sub>o\<^sub>p)"
+    using op_norm_matrix_matrix_mult_le by blast
+  finally have "((\<parallel>A\<parallel>\<^sub>o\<^sub>p))\<^sup>2 \<le> (\<parallel>transpose A\<parallel>\<^sub>o\<^sub>p) * (\<parallel>A\<parallel>\<^sub>o\<^sub>p)"
+    by linarith
+  thus "(\<parallel>A\<parallel>\<^sub>o\<^sub>p) \<le> (\<parallel>transpose A\<parallel>\<^sub>o\<^sub>p)"
+    using sq_le_cancel[of "(\<parallel>A\<parallel>\<^sub>o\<^sub>p)"] op_norm_ge_0 by metis
+qed
+
+
+subsection\<open> Matrix maximum norm \<close>
+
+abbreviation max_norm :: "real^'n^'m \<Rightarrow> real" ("(1\<parallel>_\<parallel>\<^sub>m\<^sub>a\<^sub>x)" [65] 61)
+  where "\<parallel>A\<parallel>\<^sub>m\<^sub>a\<^sub>x \<equiv> Max (abs ` (entries A))"
+
+lemma max_norm_def: "\<parallel>A\<parallel>\<^sub>m\<^sub>a\<^sub>x = Max {\<bar>A $ i $ j\<bar>|i j. i\<in>UNIV \<and> j\<in>UNIV}"
+  by (simp add: image_def, rule arg_cong[of _ _ Max], blast)
+
+lemma max_norm_set_proptys: "finite {\<bar>A $ i $ j\<bar> |i j. i \<in> UNIV \<and> j \<in> UNIV}" (is "finite ?X")
+proof-
+  have "\<And>i. finite {\<bar>A $ i $ j\<bar> | j. j \<in> UNIV}"
+    using finite_Atleast_Atmost_nat by fastforce
+  hence "finite (\<Union>i\<in>UNIV. {\<bar>A $ i $ j\<bar> | j. j \<in> UNIV})" (is "finite ?Y")
+    using finite_class.finite_UNIV by blast
+  also have "?X \<subseteq> ?Y" 
+    by auto
+  ultimately show ?thesis 
+    using finite_subset by blast
+qed
+
+lemma max_norm_ge_0: "0 \<le> \<parallel>A\<parallel>\<^sub>m\<^sub>a\<^sub>x"
+  unfolding max_norm_def 
+  apply(rule order.trans[OF abs_ge_zero[of "A $ _ $ _"] Max_ge])
+  using max_norm_set_proptys by auto
+
+lemma op_norm_le_max_norm:
+  fixes A :: "real^('n::finite)^('m::finite)"
+  shows "\<parallel>A\<parallel>\<^sub>o\<^sub>p \<le> real CARD('m) * real CARD('n) * (\<parallel>A\<parallel>\<^sub>m\<^sub>a\<^sub>x)"
+  apply(rule onorm_le_matrix_component)
+  unfolding max_norm_def by(rule Max_ge[OF max_norm_set_proptys]) force
+
+lemma sqrt_Sup_power2_eq_Sup_abs:
+  "finite A \<Longrightarrow> A \<noteq> {} \<Longrightarrow> sqrt (Sup {(f i)\<^sup>2 |i. i \<in> A}) = Sup {\<bar>f i\<bar> |i. i \<in> A}"
+proof(rule sym)
+  assume assms: "finite A" "A \<noteq> {}"
+  then obtain i where i_def: "i \<in> A \<and> Sup {(f i)\<^sup>2|i. i \<in> A} = (f i)^2"
+    using cSup_finite_ex[of "{(f i)\<^sup>2|i. i \<in> A}"] by auto
+  hence lhs: "sqrt (Sup {(f i)\<^sup>2 |i. i \<in> A}) = \<bar>f i\<bar>"
+    by simp
+  have "finite {(f i)\<^sup>2|i. i \<in> A}"
+    using assms by simp
+  hence "\<forall>j\<in>A. (f j)\<^sup>2 \<le> (f i)\<^sup>2"
+    using i_def cSup_upper[of _ "{(f i)\<^sup>2 |i. i \<in> A}"] by force
+  hence "\<forall>j\<in>A. \<bar>f j\<bar> \<le> \<bar>f i\<bar>"
+    using abs_le_square_iff by blast
+  also have "\<bar>f i\<bar> \<in> {\<bar>f i\<bar> |i. i \<in> A}"
+    using i_def by auto
+  ultimately show "Sup {\<bar>f i\<bar> |i. i \<in> A} = sqrt (Sup {(f i)\<^sup>2 |i. i \<in> A})"
+    using cSup_mem_eq[of "\<bar>f i\<bar>" "{\<bar>f i\<bar> |i. i \<in> A}"] lhs by auto
+qed
+
+lemma sqrt_Max_power2_eq_max_abs:
+  "finite A \<Longrightarrow> A \<noteq> {} \<Longrightarrow> sqrt (Max {(f i)\<^sup>2|i. i \<in> A}) = Max {\<bar>f i\<bar> |i. i \<in> A}"
+  apply(subst cSup_eq_Max[symmetric], simp_all)+
+  using sqrt_Sup_power2_eq_Sup_abs .
+
+lemma op_norm_diag_mat_eq: "\<parallel>diag_mat f\<parallel>\<^sub>o\<^sub>p = Max {\<bar>f i\<bar> |i. i \<in> UNIV}" (is "_ = Max ?A")
+proof(unfold op_norm_def)
+  have obs: "\<And>x i. (f i)\<^sup>2 * (x $ i)\<^sup>2 \<le> Max {(f i)\<^sup>2|i. i \<in> UNIV} * (x $ i)\<^sup>2"
+    apply(rule mult_right_mono[OF _ zero_le_power2])
+    using le_max_image_of_finite[of "\<lambda>i. (f i)^2"] by simp
+  {fix r assume "r \<in> {\<parallel>diag_mat f *v x\<parallel> |x. \<parallel>x\<parallel> = 1}"
+    then obtain x where x_def: "\<parallel>diag_mat f *v x\<parallel> = r \<and> \<parallel>x\<parallel> = 1"
+      by blast
+    hence "r\<^sup>2 = (\<Sum>i\<in>UNIV. (f i)\<^sup>2 * (x $ i)\<^sup>2)"
+      unfolding norm_vec_def L2_set_def matrix_vector_mul_diag_mat 
+      apply (simp add: power_mult_distrib)
+      by (metis (no_types, lifting) x_def norm_ge_zero real_sqrt_ge_0_iff real_sqrt_pow2)
+    also have "... \<le> (Max {(f i)\<^sup>2|i. i \<in> UNIV}) * (\<Sum>i\<in>UNIV. (x $ i)\<^sup>2)"
+      using obs[of _ x] by (simp add: sum_mono sum_distrib_left)
+    also have "... = Max {(f i)\<^sup>2|i. i \<in> UNIV}"
+      using x_def by (simp add: norm_vec_def L2_set_def)
+    finally have "r \<le> sqrt (Max {(f i)\<^sup>2|i. i \<in> UNIV})"
+      using x_def real_le_rsqrt by blast 
+    hence "r \<le> Max ?A"
+      by (subst (asm) sqrt_Max_power2_eq_max_abs[of UNIV f], simp_all)}
+  hence 1: "\<forall>x\<in>{\<parallel>diag_mat f *v x\<parallel> |x. \<parallel>x\<parallel> = 1}. x \<le> Max ?A"
+    unfolding diag_mat_def by blast
+  obtain i where i_def: "Max ?A = \<parallel>diag_mat f *v \<e> i\<parallel>"
+    using cMax_finite_ex[of ?A] by force
+  hence 2: "\<exists>x\<in>{\<parallel>diag_mat f *v x\<parallel> |x. \<parallel>x\<parallel> = 1}. Max ?A \<le> x"
+    by (metis (mono_tags, lifting) abs_1 mem_Collect_eq norm_axis_eq order_refl real_norm_def)
+  show "Sup {\<parallel>diag_mat f *v x\<parallel> |x. \<parallel>x\<parallel> = 1} = Max ?A"
+    by (rule cSup_eq[OF 1 2])
+qed
+
+lemma op_max_norms_eq_at_diag: "\<parallel>diag_mat f\<parallel>\<^sub>o\<^sub>p = \<parallel>diag_mat f\<parallel>\<^sub>m\<^sub>a\<^sub>x"
+proof(rule antisym)
+  have "{\<bar>f i\<bar> |i. i \<in> UNIV} \<subseteq> {\<bar>diag_mat f $ i $ j\<bar> |i j. i \<in> UNIV \<and> j \<in> UNIV}"
+    by (smt Collect_mono diag_mat_vec_nth_simps(1))
+  thus "\<parallel>diag_mat f\<parallel>\<^sub>o\<^sub>p \<le> \<parallel>diag_mat f\<parallel>\<^sub>m\<^sub>a\<^sub>x"
+    unfolding op_norm_diag_mat_eq max_norm_def
+    by (rule Max.subset_imp) (blast, simp only: finite_image_of_finite2)
+next
+  have "Sup {\<bar>diag_mat f $ i $ j\<bar> |i j. i \<in> UNIV \<and> j \<in> UNIV} \<le> Sup {\<bar>f i\<bar> |i. i \<in> UNIV}"
+    apply(rule cSup_least, blast, clarify, case_tac "i = j", simp)
+    by (rule cSup_upper, blast, simp_all) (rule cSup_upper2, auto)
+  thus "\<parallel>diag_mat f\<parallel>\<^sub>m\<^sub>a\<^sub>x \<le> \<parallel>diag_mat f\<parallel>\<^sub>o\<^sub>p"
+    unfolding op_norm_diag_mat_eq max_norm_def
+    apply (subst cSup_eq_Max[symmetric], simp only: finite_image_of_finite2, blast)
+    by (subst cSup_eq_Max[symmetric], simp, blast)
+qed
+
+
+end
\ No newline at end of file
diff --git a/thys/Matrices_for_ODEs/MTX_Preliminaries.thy b/thys/Matrices_for_ODEs/MTX_Preliminaries.thy
new file mode 100644
--- /dev/null
+++ b/thys/Matrices_for_ODEs/MTX_Preliminaries.thy
@@ -0,0 +1,403 @@
+(*  Title:       Mathematical Preliminaries
+    Author:      Jonathan Julián Huerta y Munive, 2020
+    Maintainer:  Jonathan Julián Huerta y Munive <jjhuertaymunive1@sheffield.ac.uk>
+*)
+
+section \<open> Mathematical Preliminaries \<close>
+
+text \<open>This section adds useful syntax, abbreviations and theorems to the Isabelle distribution. \<close>
+
+theory MTX_Preliminaries
+  imports "Hybrid_Systems_VCs.HS_Preliminaries"
+
+begin
+
+
+subsection \<open> Syntax \<close>
+
+abbreviation "\<e> k \<equiv> axis k 1"
+
+syntax
+  "_ivl_integral" :: "real \<Rightarrow> real \<Rightarrow> 'a \<Rightarrow> pttrn \<Rightarrow> bool" ("(3\<integral>\<^sub>_\<^sup>_ (_)\<partial>/_)" [0, 0, 10] 10)
+
+translations
+  "\<integral>\<^sub>a\<^sup>b f \<partial>x" \<rightleftharpoons> "CONST ivl_integral a b (\<lambda>x. f)"
+
+notation matrix_inv ("_\<^sup>-\<^sup>1" [90])
+
+abbreviation "entries (A::'a^'n^'m) \<equiv> {A $ i $ j | i j. i \<in> UNIV \<and> j \<in> UNIV}"
+
+
+subsection \<open> Topology and sets \<close>
+
+lemmas compact_imp_bdd_above = compact_imp_bounded[THEN bounded_imp_bdd_above]
+
+lemma comp_cont_image_spec: "continuous_on T f \<Longrightarrow> compact T \<Longrightarrow> compact {f t |t. t \<in> T}"
+  using compact_continuous_image by (simp add: Setcompr_eq_image)
+
+lemmas bdd_above_cont_comp_spec = compact_imp_bdd_above[OF comp_cont_image_spec]
+
+lemmas bdd_above_norm_cont_comp = continuous_on_norm[THEN bdd_above_cont_comp_spec]
+
+lemma open_cballE: "t\<^sub>0 \<in> T \<Longrightarrow> open T \<Longrightarrow> \<exists>e>0. cball t\<^sub>0 e \<subseteq> T"
+  using open_contains_cball by blast
+
+lemma open_ballE: "t\<^sub>0 \<in> T \<Longrightarrow> open T \<Longrightarrow> \<exists>e>0. ball t\<^sub>0 e \<subseteq> T"
+  using open_contains_ball by blast
+
+lemma funcset_UNIV: "f \<in> A \<rightarrow> UNIV"
+  by auto
+
+lemma finite_image_of_finite[simp]:
+  fixes f::"'a::finite \<Rightarrow> 'b"
+  shows "finite {x. \<exists>i. x = f i}"
+  using finite_Atleast_Atmost_nat by force
+
+lemma finite_image_of_finite2:
+  fixes f :: "'a::finite \<Rightarrow> 'b::finite \<Rightarrow> 'c"
+  shows "finite {f x y |x y. P x y}"
+proof-
+  have "finite (\<Union>x. {f x y|y. P x y})"
+    by simp
+  moreover have "{f x y|x y. P x y} = (\<Union>x. {f x y|y. P x y})"
+    by auto
+  ultimately show ?thesis 
+    by simp
+qed
+
+
+subsection \<open> Functions \<close>
+
+lemma finite_sum_univ_singleton: "(sum g UNIV) = sum g {i::'a::finite} + sum g (UNIV - {i})"
+  by (metis add.commute finite_class.finite_UNIV sum.subset_diff top_greatest)
+
+lemma suminfI:
+  fixes f :: "nat \<Rightarrow> 'a::{t2_space,comm_monoid_add}"
+  shows "f sums k \<Longrightarrow> suminf f = k"
+  unfolding sums_iff by simp
+
+lemma suminf_eq_sum:
+  fixes f :: "nat \<Rightarrow> ('a::real_normed_vector)"
+  assumes "\<And>n. n > m \<Longrightarrow> f n = 0"
+  shows "(\<Sum>n. f n) = (\<Sum>n \<le> m. f n)"
+  using assms by (meson atMost_iff finite_atMost not_le suminf_finite)
+
+lemma suminf_multr: "summable f \<Longrightarrow> (\<Sum>n. f n * c) = (\<Sum>n. f n) * c" for c::"'a::real_normed_algebra"
+  by (rule bounded_linear.suminf [OF bounded_linear_mult_left, symmetric])
+
+lemma sum_if_then_else_simps[simp]:
+  fixes q :: "('a::semiring_0)" and i :: "'n::finite"
+  shows "(\<Sum>j\<in>UNIV. f j * (if j = i then q else 0)) = f i * q"
+    and "(\<Sum>j\<in>UNIV. f j * (if i = j then q else 0)) = f i * q"
+    and "(\<Sum>j\<in>UNIV. (if i = j then q else 0) * f j) = q * f i"
+    and "(\<Sum>j\<in>UNIV. (if j = i then q else 0) * f j) = q * f i"
+  by (auto simp: finite_sum_univ_singleton[of _ i])
+
+
+subsection \<open> Suprema \<close>
+
+lemma le_max_image_of_finite[simp]:
+  fixes f::"'a::finite \<Rightarrow> 'b::linorder"
+  shows "(f i) \<le> Max {x. \<exists>i. x = f i}"
+  by (rule Max.coboundedI, simp_all) (rule_tac x=i in exI, simp)
+
+lemma cSup_eq:
+  fixes c::"'a::conditionally_complete_lattice"
+  assumes "\<forall>x \<in> X. x \<le> c" and "\<exists>x \<in> X. c \<le> x"
+  shows "Sup X = c"
+  by (metis assms cSup_eq_maximum order_class.order.antisym)
+
+lemma cSup_mem_eq: 
+  "c \<in> X \<Longrightarrow> \<forall>x \<in> X. x \<le> c \<Longrightarrow> Sup X = c" for c::"'a::conditionally_complete_lattice"
+  by (rule cSup_eq, auto)
+
+lemma cSup_finite_ex:
+  "finite X \<Longrightarrow> X \<noteq> {} \<Longrightarrow> \<exists>x\<in>X. Sup X = x" for X::"'a::conditionally_complete_linorder set"
+  by (metis (full_types) bdd_finite(1) cSup_upper finite_Sup_less_iff order_less_le)
+
+lemma cMax_finite_ex:
+  "finite X \<Longrightarrow> X \<noteq> {} \<Longrightarrow> \<exists>x\<in>X. Max X = x" for X::"'a::conditionally_complete_linorder set"
+  apply(subst cSup_eq_Max[symmetric])
+  using cSup_finite_ex by auto
+
+lemma finite_nat_minimal_witness:
+  fixes P :: "('a::finite) \<Rightarrow> nat \<Rightarrow> bool"
+  assumes "\<forall>i. \<exists>N::nat. \<forall>n \<ge> N. P i n"
+  shows "\<exists>N. \<forall>i. \<forall>n \<ge> N. P i n" 
+proof-
+  let "?bound i" = "(LEAST N. \<forall>n \<ge> N. P i n)"
+  let ?N = "Max {?bound i |i. i \<in> UNIV}"
+  {fix n::nat and i::'a 
+    assume "n \<ge> ?N" 
+    obtain M where "\<forall>n \<ge> M. P i n" 
+      using assms by blast
+    hence obs: "\<forall> m \<ge> ?bound i. P i m"
+      using LeastI[of "\<lambda>N. \<forall>n \<ge> N. P i n"] by blast
+    have "finite {?bound i |i. i \<in> UNIV}"
+      by simp
+    hence "?N \<ge> ?bound i"
+      using Max_ge by blast
+    hence "n \<ge> ?bound i" 
+      using \<open>n \<ge> ?N\<close> by linarith
+    hence "P i n" 
+      using obs by blast}
+  thus "\<exists>N. \<forall>i n. N \<le> n \<longrightarrow> P i n" 
+    by blast
+qed
+
+
+subsection \<open> Real numbers \<close>
+
+named_theorems field_power_simps "simplification rules for powers to the nth"
+
+declare semiring_normalization_rules(18) [field_power_simps]
+    and semiring_normalization_rules(26) [field_power_simps]
+    and semiring_normalization_rules(27) [field_power_simps]
+    and semiring_normalization_rules(28) [field_power_simps]
+    and semiring_normalization_rules(29) [field_power_simps]
+
+text \<open>WARNING: Adding @{thm semiring_normalization_rules(27)} to our tactic makes 
+its combination with simp to loop infinitely in some proofs.\<close>
+
+lemma sq_le_cancel:
+  shows "(a::real) \<ge> 0 \<Longrightarrow> b \<ge> 0 \<Longrightarrow> a^2 \<le> b * a \<Longrightarrow> a \<le> b"
+  and "(a::real) \<ge> 0 \<Longrightarrow> b \<ge> 0 \<Longrightarrow> a^2 \<le> a * b \<Longrightarrow> a \<le> b"
+   apply (metis less_eq_real_def mult.commute mult_le_cancel_left semiring_normalization_rules(29))
+  by (metis less_eq_real_def mult_le_cancel_left semiring_normalization_rules(29))
+
+lemma frac_diff_eq1: "a \<noteq> b \<Longrightarrow> a / (a - b) - b / (a - b) = 1" for a::real
+  by (metis (no_types, hide_lams) ab_left_minus add.commute add_left_cancel 
+      diff_divide_distrib diff_minus_eq_add div_self)
+
+lemma exp_add: "x * y - y * x = 0 \<Longrightarrow> exp (x + y) = exp x * exp y" 
+  by (rule exp_add_commuting) (simp add: ac_simps)
+
+lemmas mult_exp_exp = exp_add[symmetric]
+
+
+subsection \<open> Vectors and matrices \<close>
+
+lemma sum_axis[simp]:
+  fixes q :: "('a::semiring_0)"
+  shows "(\<Sum>j\<in>UNIV. f j * axis i q $ j) = f i * q"
+    and "(\<Sum>j\<in>UNIV. axis i q $ j * f j) = q * f i"
+  unfolding axis_def by(auto simp: vec_eq_iff)
+
+lemma sum_scalar_nth_axis: "sum (\<lambda>i. (x $ i) *s \<e> i) UNIV = x" for x :: "('a::semiring_1)^'n"
+  unfolding vec_eq_iff axis_def by simp
+
+lemma scalar_eq_scaleR[simp]: "c *s x = c *\<^sub>R x"
+  unfolding vec_eq_iff by simp
+
+lemma matrix_add_rdistrib: "((B + C) ** A) = (B ** A) + (C ** A)"
+  by (vector matrix_matrix_mult_def sum.distrib[symmetric] field_simps)
+
+lemma vec_mult_inner: "(A *v v) \<bullet> w = v \<bullet> (transpose A *v w)" for A :: "real ^'n^'n"
+  unfolding matrix_vector_mult_def transpose_def inner_vec_def
+  apply(simp add: sum_distrib_right sum_distrib_left)
+  apply(subst sum.swap)
+  apply(subgoal_tac "\<forall>i j. A $ i $ j * v $ j * w $ i = v $ j * (A $ i $ j * w $ i)")
+  by presburger simp
+
+lemma uminus_axis_eq[simp]: "- axis i k = axis i (-k)" for k :: "'a::ring"
+  unfolding axis_def by(simp add: vec_eq_iff)
+
+lemma norm_axis_eq[simp]: "\<parallel>axis i k\<parallel> = \<parallel>k\<parallel>"
+proof(simp add: axis_def norm_vec_def L2_set_def)
+  let "?\<delta>\<^sub>K" = "\<lambda>i j k. if i = j then k else 0" 
+  have "(\<Sum>j\<in>UNIV. (\<parallel>(?\<delta>\<^sub>K j i k)\<parallel>)\<^sup>2) = (\<Sum>j\<in>{i}. (\<parallel>(?\<delta>\<^sub>K j i k)\<parallel>)\<^sup>2) + (\<Sum>j\<in>(UNIV-{i}). (\<parallel>(?\<delta>\<^sub>K j i k)\<parallel>)\<^sup>2)"
+    using finite_sum_univ_singleton by blast 
+  also have "... = (\<parallel>k\<parallel>)\<^sup>2" 
+    by simp
+  finally show "sqrt (\<Sum>j\<in>UNIV. (norm (if j = i then k else 0))\<^sup>2) = norm k" 
+    by simp
+qed
+
+lemma matrix_axis_0: 
+  fixes A :: "('a::idom)^'n^'m"
+  assumes "k \<noteq> 0 " and h:"\<forall>i. (A *v (axis i k)) = 0"
+  shows "A = 0"  
+proof-
+  {fix i::'n
+    have "0 = (\<Sum>j\<in>UNIV. (axis i k) $ j *s column j A)" 
+      using h matrix_mult_sum[of A "axis i k"] by simp
+    also have "... = k *s column i A" 
+      by (simp add: axis_def vector_scalar_mult_def column_def vec_eq_iff mult.commute)
+    finally have "k *s column i A = 0"
+      unfolding axis_def by simp
+    hence "column i A = 0" 
+      using vector_mul_eq_0 \<open>k \<noteq> 0\<close> by blast}
+  thus "A = 0" 
+    unfolding column_def vec_eq_iff by simp
+qed
+
+lemma scaleR_norm_sgn_eq: "(\<parallel>x\<parallel>) *\<^sub>R sgn x = x"
+  by (metis divideR_right norm_eq_zero scale_eq_0_iff sgn_div_norm)
+
+lemma vector_scaleR_commute: "A *v c *\<^sub>R x = c *\<^sub>R (A *v x)" for x :: "('a::real_normed_algebra_1)^'n"
+  unfolding scaleR_vec_def matrix_vector_mult_def by(auto simp: vec_eq_iff scaleR_right.sum)
+
+lemma scaleR_vector_assoc: "c *\<^sub>R (A *v x) = (c *\<^sub>R A) *v x" for x :: "('a::real_normed_algebra_1)^'n"
+  unfolding matrix_vector_mult_def by(auto simp: vec_eq_iff scaleR_right.sum)
+
+lemma mult_norm_matrix_sgn_eq:
+  fixes x :: "('a::real_normed_algebra_1)^'n"
+  shows "(\<parallel>A *v sgn x\<parallel>) * (\<parallel>x\<parallel>) = \<parallel>A *v x\<parallel>"
+proof-
+  have "\<parallel>A *v x\<parallel> = \<parallel>A *v ((\<parallel>x\<parallel>) *\<^sub>R sgn x)\<parallel>"
+    by(simp add: scaleR_norm_sgn_eq)
+  also have "... = (\<parallel>A *v sgn x\<parallel>) * (\<parallel>x\<parallel>)"
+    by(simp add: vector_scaleR_commute)
+  finally show ?thesis ..
+qed
+
+
+subsection\<open> Diagonalization \<close>
+
+lemma invertibleI: "A ** B = mat 1 \<Longrightarrow> B ** A = mat 1 \<Longrightarrow> invertible A"
+  unfolding invertible_def by auto
+
+lemma invertibleD[simp]:
+  assumes "invertible A" 
+  shows "A\<^sup>-\<^sup>1 ** A = mat 1" and "A ** A\<^sup>-\<^sup>1 = mat 1"
+  using assms unfolding matrix_inv_def invertible_def
+  by (simp_all add: verit_sko_ex')
+
+lemma matrix_inv_unique:
+  assumes "A ** B = mat 1" and "B ** A = mat 1"
+  shows "A\<^sup>-\<^sup>1 = B"
+  by (metis assms invertibleD(2) invertibleI matrix_mul_assoc matrix_mul_lid)
+
+lemma invertible_matrix_inv: "invertible A \<Longrightarrow> invertible (A\<^sup>-\<^sup>1)"
+  using invertibleD invertibleI by blast
+
+lemma matrix_inv_idempotent[simp]: "invertible A \<Longrightarrow> A\<^sup>-\<^sup>1\<^sup>-\<^sup>1 = A"
+  using invertibleD matrix_inv_unique by blast
+  
+lemma matrix_inv_matrix_mul:
+  assumes "invertible A" and "invertible B"
+  shows "(A ** B)\<^sup>-\<^sup>1 = B\<^sup>-\<^sup>1 ** A\<^sup>-\<^sup>1"
+proof(rule matrix_inv_unique)
+  have "A ** B ** (B\<^sup>-\<^sup>1 ** A\<^sup>-\<^sup>1) = A ** (B ** B\<^sup>-\<^sup>1) ** A\<^sup>-\<^sup>1"
+    by (simp add: matrix_mul_assoc)
+  also have "... = mat 1"
+    using assms by simp
+  finally show "A ** B ** (B\<^sup>-\<^sup>1 ** A\<^sup>-\<^sup>1) = mat 1" .
+next
+  have "B\<^sup>-\<^sup>1 ** A\<^sup>-\<^sup>1 ** (A ** B) = B\<^sup>-\<^sup>1 ** (A\<^sup>-\<^sup>1 ** A) ** B"
+    by (simp add: matrix_mul_assoc)
+  also have "... = mat 1"
+    using assms by simp
+  finally show "B\<^sup>-\<^sup>1 ** A\<^sup>-\<^sup>1 ** (A ** B) = mat 1" .
+qed
+
+lemma mat_inverse_simps[simp]:
+  fixes c :: "'a::division_ring"
+  assumes "c \<noteq> 0"
+  shows "mat (inverse c) ** mat c = mat 1" 
+    and "mat c ** mat (inverse c) = mat 1"
+  unfolding matrix_matrix_mult_def mat_def by (auto simp: vec_eq_iff assms)
+
+lemma matrix_inv_mat[simp]: "c \<noteq> 0 \<Longrightarrow> (mat c)\<^sup>-\<^sup>1 = mat (inverse c)" for c :: "'a::division_ring"
+  by (simp add: matrix_inv_unique)
+
+lemma invertible_mat[simp]: "c \<noteq> 0 \<Longrightarrow> invertible (mat c)" for c :: "'a::division_ring"
+  using invertibleI mat_inverse_simps(1) mat_inverse_simps(2) by blast
+
+lemma matrix_inv_mat_1: "(mat (1::'a::division_ring))\<^sup>-\<^sup>1 = mat 1"
+  by simp
+
+lemma invertible_mat_1: "invertible (mat (1::'a::division_ring))"
+  by simp
+
+definition similar_matrix :: "('a::semiring_1)^'m^'m \<Rightarrow> ('a::semiring_1)^'n^'n \<Rightarrow> bool" (infixr "\<sim>" 25)
+  where "similar_matrix A B \<longleftrightarrow> (\<exists> P. invertible P \<and> A = P\<^sup>-\<^sup>1 ** B ** P)"
+
+lemma similar_matrix_refl[simp]: "A \<sim> A" for A :: "'a::division_ring^'n^'n"
+  by (unfold similar_matrix_def, rule_tac x="mat 1" in exI, simp)
+
+lemma similar_matrix_simm: "A \<sim> B \<Longrightarrow> B \<sim> A" for A B :: "('a::semiring_1)^'n^'n"
+  apply(unfold similar_matrix_def, clarsimp)
+  apply(rule_tac x="P\<^sup>-\<^sup>1" in exI, simp add: invertible_matrix_inv)
+  by (metis invertible_def matrix_inv_unique matrix_mul_assoc matrix_mul_lid matrix_mul_rid)
+
+lemma similar_matrix_trans: "A \<sim> B \<Longrightarrow> B \<sim> C \<Longrightarrow> A \<sim> C" for A B C :: "('a::semiring_1)^'n^'n"
+proof(unfold similar_matrix_def, clarsimp)
+  fix P Q
+  assume "A = P\<^sup>-\<^sup>1 ** (Q\<^sup>-\<^sup>1 ** C ** Q) ** P" and "B = Q\<^sup>-\<^sup>1 ** C ** Q"
+  let ?R = "Q ** P"
+  assume inverts: "invertible Q" "invertible P"
+  hence "?R\<^sup>-\<^sup>1 = P\<^sup>-\<^sup>1 ** Q\<^sup>-\<^sup>1"
+    by (rule matrix_inv_matrix_mul)
+  also have "invertible ?R"
+    using inverts invertible_mult by blast 
+  ultimately show "\<exists>R. invertible R \<and> P\<^sup>-\<^sup>1 ** (Q\<^sup>-\<^sup>1 ** C ** Q) ** P = R\<^sup>-\<^sup>1 ** C ** R"
+    by (metis matrix_mul_assoc)
+qed
+
+lemma mat_vec_nth_simps[simp]:
+  "i = j \<Longrightarrow> mat c $ i $ j = c" 
+  "i \<noteq> j \<Longrightarrow> mat c $ i $ j = 0"
+  by (simp_all add: mat_def)
+
+definition "diag_mat f = (\<chi> i j. if i = j then f i else 0)"
+
+lemma diag_mat_vec_nth_simps[simp]:
+  "i = j \<Longrightarrow> diag_mat f $ i $ j = f i"
+  "i \<noteq> j \<Longrightarrow> diag_mat f $ i $ j = 0"
+  unfolding diag_mat_def by simp_all
+
+lemma diag_mat_const_eq[simp]: "diag_mat (\<lambda>i. c) = mat c"
+  unfolding mat_def diag_mat_def by simp
+
+lemma matrix_vector_mul_diag_mat: "diag_mat f *v s = (\<chi> i. f i * s$i)"
+  unfolding diag_mat_def matrix_vector_mult_def by simp
+
+lemma matrix_vector_mul_diag_axis[simp]: "diag_mat f *v (axis i k) = axis i (f i * k)"
+  by (simp add: matrix_vector_mul_diag_mat axis_def fun_eq_iff)
+
+lemma matrix_mul_diag_matl: "diag_mat f ** A = (\<chi> i j. f i * A$i$j)"
+  unfolding diag_mat_def matrix_matrix_mult_def by simp
+
+lemma matrix_matrix_mul_diag_matr: "A ** diag_mat f = (\<chi> i j. A$i$j * f j)"
+  unfolding diag_mat_def matrix_matrix_mult_def apply(clarsimp simp: fun_eq_iff)
+  subgoal for i j 
+    by (auto simp: finite_sum_univ_singleton[of _ j])
+  done
+
+lemma matrix_mul_diag_diag: "diag_mat f ** diag_mat g = diag_mat (\<lambda>i. f i * g i)"
+  unfolding diag_mat_def matrix_matrix_mult_def vec_eq_iff by simp
+
+lemma compow_matrix_mul_diag_mat_eq: "((**) (diag_mat f) ^^ n) (mat 1) = diag_mat (\<lambda>i. f i^n)"
+  apply(induct n, simp_all add: matrix_mul_diag_matl)
+  by (auto simp: vec_eq_iff diag_mat_def)
+
+lemma compow_similar_diag_mat_eq:
+  assumes "invertible P" 
+      and "A = P\<^sup>-\<^sup>1 ** (diag_mat f) ** P"
+    shows "((**) A ^^ n) (mat 1) = P\<^sup>-\<^sup>1 ** (diag_mat (\<lambda>i. f i^n)) ** P"
+proof(induct n, simp_all add: assms)
+  fix n::nat
+  have "P\<^sup>-\<^sup>1 ** diag_mat f ** P ** (P\<^sup>-\<^sup>1 ** diag_mat (\<lambda>i. f i ^ n) ** P) = 
+  P\<^sup>-\<^sup>1 ** diag_mat f ** diag_mat (\<lambda>i. f i ^ n) ** P" (is "?lhs = _")
+    by (metis (no_types, lifting) assms(1) invertibleD(2) matrix_mul_rid matrix_mul_assoc) 
+  also have "... = P\<^sup>-\<^sup>1 ** diag_mat (\<lambda>i. f i * f i ^ n) ** P" (is "_ = ?rhs")
+    by (metis (full_types) matrix_mul_assoc matrix_mul_diag_diag)
+  finally show "?lhs = ?rhs" .
+qed
+
+lemma compow_similar_diag_mat:
+  assumes "A \<sim> (diag_mat f)"
+  shows "((**) A ^^ n) (mat 1) \<sim> diag_mat (\<lambda>i. f i^n)"
+proof(unfold similar_matrix_def)
+  obtain P where "invertible P" and "A = P\<^sup>-\<^sup>1 ** (diag_mat f) ** P"
+    using assms unfolding similar_matrix_def by blast
+  thus "\<exists>P. invertible P \<and> ((**) A ^^ n) (mat 1) = P\<^sup>-\<^sup>1 ** diag_mat (\<lambda>i. f i ^ n) ** P"
+    using compow_similar_diag_mat_eq by blast
+qed
+
+no_notation matrix_inv ("_\<^sup>-\<^sup>1" [90])
+        and similar_matrix (infixr "\<sim>" 25)
+
+
+end
\ No newline at end of file
diff --git a/thys/Matrices_for_ODEs/ROOT b/thys/Matrices_for_ODEs/ROOT
new file mode 100644
--- /dev/null
+++ b/thys/Matrices_for_ODEs/ROOT
@@ -0,0 +1,11 @@
+chapter AFP
+
+session "Matrices_for_ODEs" (AFP) = "HOL-Analysis" +
+  options [timeout = 1800]
+  sessions
+    Hybrid_Systems_VCs
+  theories
+    MTX_Examples
+  document_files
+    "root.bib"
+    "root.tex"
\ No newline at end of file
diff --git a/thys/Matrices_for_ODEs/SQ_MTX.thy b/thys/Matrices_for_ODEs/SQ_MTX.thy
new file mode 100644
--- /dev/null
+++ b/thys/Matrices_for_ODEs/SQ_MTX.thy
@@ -0,0 +1,718 @@
+(*  Title:       Square Matrices
+    Author:      Jonathan Julián Huerta y Munive, 2020
+    Maintainer:  Jonathan Julián Huerta y Munive <jjhuertaymunive1@sheffield.ac.uk>
+*)
+
+section \<open> Square Matrices \<close>
+
+text\<open> The general solution for affine systems of ODEs involves the exponential function. 
+Unfortunately, this operation is only available in Isabelle for the type class ``banach''. 
+Hence, we define a type of square matrices and prove that it is an instance of this class.\<close>
+
+theory SQ_MTX
+  imports MTX_Norms
+
+begin
+
+subsection \<open> Definition \<close>
+
+typedef 'm sq_mtx = "UNIV::(real^'m^'m) set"
+  morphisms to_vec to_mtx by simp
+
+declare to_mtx_inverse [simp]
+    and to_vec_inverse [simp]
+
+setup_lifting type_definition_sq_mtx
+
+lift_definition sq_mtx_ith :: "'m sq_mtx \<Rightarrow> 'm \<Rightarrow> (real^'m)" (infixl "$$" 90) is "($)" .
+
+lift_definition sq_mtx_vec_mult :: "'m sq_mtx \<Rightarrow> (real^'m) \<Rightarrow> (real^'m)" (infixl "*\<^sub>V" 90) is "(*v)" .
+
+lift_definition vec_sq_mtx_prod :: "(real^'m) \<Rightarrow> 'm sq_mtx \<Rightarrow> (real^'m)" is "(v*)" .
+
+lift_definition sq_mtx_diag :: "(('m::finite) \<Rightarrow> real) \<Rightarrow> ('m::finite) sq_mtx" (binder "\<d>\<i>\<a>\<g> " 10) 
+  is diag_mat .
+
+lift_definition sq_mtx_transpose :: "('m::finite) sq_mtx \<Rightarrow> 'm sq_mtx" ("_\<^sup>\<dagger>") is transpose .
+
+lift_definition sq_mtx_inv :: "('m::finite) sq_mtx \<Rightarrow> 'm sq_mtx" ("_\<^sup>-\<^sup>1" [90]) is matrix_inv .
+
+lift_definition sq_mtx_row :: "'m \<Rightarrow> ('m::finite) sq_mtx \<Rightarrow> real^'m" ("\<r>\<o>\<w>") is row .
+
+lift_definition sq_mtx_col :: "'m \<Rightarrow> ('m::finite) sq_mtx \<Rightarrow> real^'m" ("\<c>\<o>\<l>")  is column .
+
+lemma to_vec_eq_ith: "(to_vec A) $ i = A $$ i"
+  by transfer simp
+
+lemma to_mtx_ith[simp]: 
+  "(to_mtx A) $$ i1 = A $ i1"
+  "(to_mtx A) $$ i1 $ i2 = A $ i1 $ i2"
+  by (transfer, simp)+
+
+lemma to_mtx_vec_lambda_ith[simp]: "to_mtx (\<chi> i j. x i j) $$ i1 $ i2 = x i1 i2"
+  by (simp add: sq_mtx_ith_def)
+
+lemma sq_mtx_eq_iff:
+  shows "A = B = (\<forall>i j. A $$ i $ j = B $$ i $ j)"
+    and "A = B = (\<forall>i. A $$ i = B $$ i)"
+  by (transfer, simp add: vec_eq_iff)+
+
+lemma sq_mtx_diag_simps[simp]:
+  "i = j \<Longrightarrow> sq_mtx_diag f $$ i $ j = f i"
+  "i \<noteq> j \<Longrightarrow> sq_mtx_diag f $$ i $ j = 0"
+  "sq_mtx_diag f $$ i = axis i (f i)"
+  unfolding sq_mtx_diag_def by (simp_all add: axis_def vec_eq_iff)
+
+lemma sq_mtx_diag_vec_mult: "(\<d>\<i>\<a>\<g> i. f i) *\<^sub>V s = (\<chi> i. f i * s$i)"
+  by (simp add: matrix_vector_mul_diag_mat sq_mtx_diag.abs_eq sq_mtx_vec_mult.abs_eq)
+
+lemma sq_mtx_vec_mult_diag_axis: "(\<d>\<i>\<a>\<g> i. f i) *\<^sub>V (axis i k) = axis i (f i * k)"
+  unfolding sq_mtx_diag_vec_mult axis_def by auto
+
+lemma sq_mtx_vec_mult_eq: "m *\<^sub>V x = (\<chi> i. sum (\<lambda>j. (m $$ i $ j) * (x $ j)) UNIV)"
+  by (transfer, simp add: matrix_vector_mult_def)
+
+lemma sq_mtx_transpose_transpose[simp]: "(A\<^sup>\<dagger>)\<^sup>\<dagger> = A"
+  by (transfer, simp)
+
+lemma transpose_mult_vec_canon_row[simp]: "(A\<^sup>\<dagger>) *\<^sub>V (\<e> i) = \<r>\<o>\<w> i A"
+  by transfer (simp add: row_def transpose_def axis_def matrix_vector_mult_def)
+
+lemma row_ith[simp]: "\<r>\<o>\<w> i A = A $$ i"
+  by transfer (simp add: row_def)
+
+lemma mtx_vec_mult_canon: "A *\<^sub>V (\<e> i) = \<c>\<o>\<l> i A" 
+  by (transfer, simp add: matrix_vector_mult_basis)
+
+
+subsection \<open> Ring of square matrices \<close>
+
+instantiation sq_mtx :: (finite) ring 
+begin
+
+lift_definition plus_sq_mtx :: "'a sq_mtx \<Rightarrow> 'a sq_mtx \<Rightarrow> 'a sq_mtx" is "(+)" .
+
+lift_definition zero_sq_mtx :: "'a sq_mtx" is "0" .
+
+lift_definition uminus_sq_mtx :: "'a sq_mtx \<Rightarrow> 'a sq_mtx" is "uminus" .
+
+lift_definition minus_sq_mtx :: "'a sq_mtx \<Rightarrow> 'a sq_mtx \<Rightarrow> 'a sq_mtx" is "(-)" .
+
+lift_definition times_sq_mtx :: "'a sq_mtx \<Rightarrow> 'a sq_mtx \<Rightarrow> 'a sq_mtx" is "(**)" .
+
+declare plus_sq_mtx.rep_eq [simp]
+    and minus_sq_mtx.rep_eq [simp]
+
+instance apply intro_classes
+  by(transfer, simp add: algebra_simps matrix_mul_assoc matrix_add_rdistrib matrix_add_ldistrib)+
+
+end
+
+lemma sq_mtx_zero_ith[simp]: "0 $$ i = 0"
+  by (transfer, simp)
+
+lemma sq_mtx_zero_nth[simp]: "0 $$ i $ j = 0"
+  by transfer simp
+
+lemma sq_mtx_plus_eq: "A + B = to_mtx (\<chi> i j. A$$i$j + B$$i$j)"
+  by transfer (simp add: vec_eq_iff)
+
+lemma sq_mtx_plus_ith[simp]:"(A + B) $$ i = A $$ i + B $$ i"
+  unfolding sq_mtx_plus_eq by (simp add: vec_eq_iff)
+
+lemma sq_mtx_uminus_eq: "- A = to_mtx (\<chi> i j. - A$$i$j)"
+  by transfer (simp add: vec_eq_iff)
+
+lemma sq_mtx_minus_eq: "A - B = to_mtx (\<chi> i j. A$$i$j - B$$i$j)"
+  by transfer (simp add: vec_eq_iff)
+
+lemma sq_mtx_minus_ith[simp]:"(A - B) $$ i = A $$ i - B $$ i"
+  unfolding sq_mtx_minus_eq by (simp add: vec_eq_iff)
+
+lemma sq_mtx_times_eq: "A * B = to_mtx (\<chi> i j. sum (\<lambda>k. A$$i$k * B$$k$j) UNIV)"
+  by transfer (simp add: matrix_matrix_mult_def)
+
+lemma sq_mtx_plus_diag_diag[simp]: "sq_mtx_diag f + sq_mtx_diag g = (\<d>\<i>\<a>\<g> i. f i + g i)"
+  by (subst sq_mtx_eq_iff) (simp add: axis_def)
+
+lemma sq_mtx_minus_diag_diag[simp]: "sq_mtx_diag f - sq_mtx_diag g = (\<d>\<i>\<a>\<g> i. f i - g i)"
+  by (subst sq_mtx_eq_iff) (simp add: axis_def)
+
+lemma sum_sq_mtx_diag[simp]: "(\<Sum>n<m. sq_mtx_diag (g n)) = (\<d>\<i>\<a>\<g> i. \<Sum>n<m. (g n i))" for m::nat
+  by (induct m, simp, subst sq_mtx_eq_iff, simp_all)
+
+lemma sq_mtx_mult_diag_diag[simp]: "sq_mtx_diag f * sq_mtx_diag g = (\<d>\<i>\<a>\<g> i. f i * g i)"
+  by (simp add: matrix_mul_diag_diag sq_mtx_diag.abs_eq times_sq_mtx.abs_eq)
+
+lemma sq_mtx_mult_diagl: "(\<d>\<i>\<a>\<g> i. f i) * A = to_mtx (\<chi> i j. f i * A $$ i $ j)"
+  by transfer (simp add: matrix_mul_diag_matl)
+
+lemma sq_mtx_mult_diagr: "A * (\<d>\<i>\<a>\<g> i. f i) = to_mtx (\<chi> i j. A $$ i $ j * f j)"
+  by transfer (simp add: matrix_matrix_mul_diag_matr)
+
+lemma mtx_vec_mult_0l[simp]: "0 *\<^sub>V x = 0"
+  by (simp add: sq_mtx_vec_mult.abs_eq zero_sq_mtx_def)
+
+lemma mtx_vec_mult_0r[simp]: "A *\<^sub>V 0 = 0"
+  by (transfer, simp)
+
+lemma mtx_vec_mult_add_rdistr: "(A + B) *\<^sub>V x = A *\<^sub>V x + B *\<^sub>V x"
+  unfolding plus_sq_mtx_def 
+  apply(transfer)
+  by (simp add: matrix_vector_mult_add_rdistrib)
+
+lemma mtx_vec_mult_add_rdistl: "A *\<^sub>V (x + y) = A *\<^sub>V x + A *\<^sub>V y"
+  unfolding plus_sq_mtx_def 
+  apply transfer
+  by (simp add: matrix_vector_right_distrib)
+
+lemma mtx_vec_mult_minus_rdistrib: "(A - B) *\<^sub>V x = A *\<^sub>V x - B *\<^sub>V x"
+  unfolding minus_sq_mtx_def by(transfer, simp add: matrix_vector_mult_diff_rdistrib)
+
+lemma mtx_vec_mult_minus_ldistrib: "A *\<^sub>V (x - y) =  A *\<^sub>V x -  A *\<^sub>V y"
+  by (metis (no_types, lifting) add_diff_cancel diff_add_cancel 
+      matrix_vector_right_distrib sq_mtx_vec_mult.rep_eq)
+
+lemma sq_mtx_times_vec_assoc: "(A * B) *\<^sub>V x = A *\<^sub>V (B *\<^sub>V x)"
+  by (transfer, simp add: matrix_vector_mul_assoc)
+
+lemma sq_mtx_vec_mult_sum_cols: "A *\<^sub>V x = sum (\<lambda>i. x $ i *\<^sub>R \<c>\<o>\<l> i A) UNIV"
+  by(transfer) (simp add: matrix_mult_sum scalar_mult_eq_scaleR)
+
+
+subsection \<open> Real normed vector space of square matrices \<close>
+
+instantiation sq_mtx :: (finite) real_normed_vector 
+begin
+
+definition norm_sq_mtx :: "'a sq_mtx \<Rightarrow> real" where "\<parallel>A\<parallel> = \<parallel>to_vec A\<parallel>\<^sub>o\<^sub>p"
+
+lift_definition scaleR_sq_mtx :: "real \<Rightarrow> 'a sq_mtx \<Rightarrow> 'a sq_mtx" is scaleR .
+
+definition sgn_sq_mtx :: "'a sq_mtx \<Rightarrow> 'a sq_mtx" 
+  where "sgn_sq_mtx A = (inverse (\<parallel>A\<parallel>)) *\<^sub>R A"
+
+definition dist_sq_mtx :: "'a sq_mtx \<Rightarrow> 'a sq_mtx \<Rightarrow> real" 
+  where "dist_sq_mtx A B = \<parallel>A - B\<parallel>" 
+
+definition uniformity_sq_mtx :: "('a sq_mtx \<times> 'a sq_mtx) filter" 
+  where "uniformity_sq_mtx = (INF e\<in>{0<..}. principal {(x, y). dist x y < e})"
+
+definition open_sq_mtx :: "'a sq_mtx set \<Rightarrow> bool" 
+  where "open_sq_mtx U = (\<forall>x\<in>U. \<forall>\<^sub>F (x', y) in uniformity. x' = x \<longrightarrow> y \<in> U)"
+
+instance apply intro_classes 
+  unfolding sgn_sq_mtx_def open_sq_mtx_def dist_sq_mtx_def uniformity_sq_mtx_def
+            prefer 10 
+            apply(transfer, simp add: norm_sq_mtx_def op_norm_triangle)
+           prefer 9 
+           apply(simp_all add: norm_sq_mtx_def zero_sq_mtx_def op_norm_eq_0)
+  by (transfer, simp add: norm_sq_mtx_def op_norm_scaleR algebra_simps)+
+
+end
+
+lemma sq_mtx_scaleR_eq: "c *\<^sub>R A = to_mtx (\<chi> i j. c *\<^sub>R A $$ i $ j)"
+  by transfer (simp add: vec_eq_iff)
+
+lemma scaleR_to_mtx_ith[simp]: "c *\<^sub>R (to_mtx A) $$ i1 $ i2 = c * A $ i1 $ i2"
+  by transfer (simp add: scaleR_vec_def)
+
+lemma sq_mtx_scaleR_ith[simp]: "(c *\<^sub>R A) $$ i = (c  *\<^sub>R (A $$ i))"
+  by (unfold scaleR_sq_mtx_def, transfer, simp)
+
+lemma scaleR_sq_mtx_diag: "c *\<^sub>R sq_mtx_diag f = (\<d>\<i>\<a>\<g> i. c * f i)"
+  by (subst sq_mtx_eq_iff, simp add: axis_def)
+
+lemma scaleR_mtx_vec_assoc: "(c *\<^sub>R A) *\<^sub>V x = c *\<^sub>R (A *\<^sub>V x)"
+  unfolding scaleR_sq_mtx_def sq_mtx_vec_mult_def apply simp
+  by (simp add: scaleR_matrix_vector_assoc)
+
+lemma mtx_vec_scaleR_commute: "A *\<^sub>V (c *\<^sub>R x) = c *\<^sub>R (A *\<^sub>V x)"
+  unfolding scaleR_sq_mtx_def sq_mtx_vec_mult_def apply(simp, transfer)
+  by (simp add: vector_scaleR_commute)
+
+lemma mtx_times_scaleR_commute: "A * (c *\<^sub>R B) = c *\<^sub>R (A * B)" for A::"('n::finite) sq_mtx"
+  unfolding sq_mtx_scaleR_eq sq_mtx_times_eq 
+  apply(simp add: to_mtx_inject)
+  apply(simp add: vec_eq_iff fun_eq_iff)
+  by (simp add: semiring_normalization_rules(19) vector_space_over_itself.scale_sum_right)
+
+lemma le_mtx_norm: "m \<in> {\<parallel>A *\<^sub>V x\<parallel> |x. \<parallel>x\<parallel> = 1} \<Longrightarrow> m \<le> \<parallel>A\<parallel>"
+  using cSup_upper[of _ "{\<parallel>(to_vec A) *v x\<parallel> | x. \<parallel>x\<parallel> = 1}"]
+  by (simp add: op_norm_set_proptys(2) op_norm_def norm_sq_mtx_def sq_mtx_vec_mult.rep_eq)
+
+lemma norm_vec_mult_le: "\<parallel>A *\<^sub>V x\<parallel> \<le> (\<parallel>A\<parallel>) * (\<parallel>x\<parallel>)"
+  by (simp add: norm_matrix_le_mult_op_norm norm_sq_mtx_def sq_mtx_vec_mult.rep_eq)
+
+lemma bounded_bilinear_sq_mtx_vec_mult: "bounded_bilinear (\<lambda>A s. A *\<^sub>V s)"
+  apply (rule bounded_bilinear.intro, simp_all add: mtx_vec_mult_add_rdistr 
+      mtx_vec_mult_add_rdistl scaleR_mtx_vec_assoc mtx_vec_scaleR_commute)
+  by (rule_tac x=1 in exI, auto intro!: norm_vec_mult_le)
+
+lemma norm_sq_mtx_def2: "\<parallel>A\<parallel> = Sup {\<parallel>A *\<^sub>V x\<parallel> |x. \<parallel>x\<parallel> = 1}"
+  unfolding norm_sq_mtx_def op_norm_def sq_mtx_vec_mult_def by simp
+
+lemma norm_sq_mtx_def3: "\<parallel>A\<parallel> = (SUP x. (\<parallel>A *\<^sub>V x\<parallel>) / (\<parallel>x\<parallel>))"
+  unfolding norm_sq_mtx_def onorm_def sq_mtx_vec_mult_def by simp
+
+lemma norm_sq_mtx_diag: "\<parallel>sq_mtx_diag f\<parallel> = Max {\<bar>f i\<bar> |i. i \<in> UNIV}"
+  unfolding norm_sq_mtx_def apply transfer
+  by (rule op_norm_diag_mat_eq)
+
+lemma sq_mtx_norm_le_sum_col: "\<parallel>A\<parallel> \<le> (\<Sum>i\<in>UNIV. \<parallel>\<c>\<o>\<l> i A\<parallel>)"
+  using op_norm_le_sum_column[of "to_vec A"] 
+  apply(simp add: norm_sq_mtx_def)
+  by(transfer, simp add: op_norm_le_sum_column)
+
+lemma norm_le_transpose: "\<parallel>A\<parallel> \<le> \<parallel>A\<^sup>\<dagger>\<parallel>"
+  unfolding norm_sq_mtx_def by transfer (rule op_norm_le_transpose)
+
+lemma norm_eq_norm_transpose[simp]: "\<parallel>A\<^sup>\<dagger>\<parallel> = \<parallel>A\<parallel>"
+  using norm_le_transpose[of A] and norm_le_transpose[of "A\<^sup>\<dagger>"] by simp
+
+lemma norm_column_le_norm: "\<parallel>A $$ i\<parallel> \<le> \<parallel>A\<parallel>"
+  using norm_vec_mult_le[of "A\<^sup>\<dagger>" "\<e> i"] by simp
+
+
+subsection \<open> Real normed algebra of square matrices \<close>
+
+instantiation sq_mtx :: (finite) real_normed_algebra_1
+begin
+
+lift_definition one_sq_mtx :: "'a sq_mtx" is "to_mtx (mat 1)" .
+
+lemma sq_mtx_one_idty: "1 * A = A" "A * 1 = A" for A :: "'a sq_mtx"
+  by(transfer, transfer, unfold mat_def matrix_matrix_mult_def, simp add: vec_eq_iff)+
+
+lemma sq_mtx_norm_1: "\<parallel>(1::'a sq_mtx)\<parallel> = 1"
+  unfolding one_sq_mtx_def norm_sq_mtx_def 
+  apply(simp add: op_norm_def)
+  apply(subst cSup_eq[of _ 1])
+  using ex_norm_eq_1 by auto
+
+lemma sq_mtx_norm_times: "\<parallel>A * B\<parallel> \<le> (\<parallel>A\<parallel>) * (\<parallel>B\<parallel>)" for A :: "'a sq_mtx"
+  unfolding norm_sq_mtx_def times_sq_mtx_def by(simp add: op_norm_matrix_matrix_mult_le)
+
+instance 
+  apply intro_classes 
+  apply(simp_all add: sq_mtx_one_idty sq_mtx_norm_1 sq_mtx_norm_times)
+  apply(simp_all add: to_mtx_inject vec_eq_iff one_sq_mtx_def zero_sq_mtx_def mat_def)
+  by(transfer, simp add: scalar_matrix_assoc matrix_scalar_ac)+
+
+end
+
+lemma sq_mtx_one_ith_simps[simp]: "1 $$ i $ i = 1" "i \<noteq> j \<Longrightarrow> 1 $$ i $ j = 0"
+  unfolding one_sq_mtx_def mat_def by simp_all
+
+lemma of_nat_eq_sq_mtx_diag[simp]: "of_nat m = (\<d>\<i>\<a>\<g> i. m)"
+  by (induct m) (simp, subst sq_mtx_eq_iff, simp add: axis_def)+
+
+lemma mtx_vec_mult_1[simp]: "1 *\<^sub>V s = s"
+  by (auto simp: sq_mtx_vec_mult_def one_sq_mtx_def 
+      mat_def vec_eq_iff matrix_vector_mult_def)
+
+lemma sq_mtx_diag_one[simp]: "(\<d>\<i>\<a>\<g> i. 1) = 1"
+  by (subst sq_mtx_eq_iff, simp add: one_sq_mtx_def mat_def axis_def)
+
+abbreviation "mtx_invertible A \<equiv> invertible (to_vec A)"
+
+lemma mtx_invertible_def: "mtx_invertible A \<longleftrightarrow> (\<exists>A'. A' * A = 1 \<and> A * A' = 1)"
+  apply (unfold sq_mtx_inv_def times_sq_mtx_def one_sq_mtx_def invertible_def, clarsimp, safe)
+   apply(rule_tac x="to_mtx A'" in exI, simp)
+  by (rule_tac x="to_vec A'" in exI, simp add: to_mtx_inject)
+
+lemma mtx_invertibleI:
+  assumes "A * B = 1" and "B * A = 1"
+  shows "mtx_invertible A"
+  using assms unfolding mtx_invertible_def by auto
+
+lemma mtx_invertibleD[simp]:
+  assumes "mtx_invertible A" 
+  shows "A\<^sup>-\<^sup>1 * A = 1" and "A * A\<^sup>-\<^sup>1 = 1"
+  apply (unfold sq_mtx_inv_def times_sq_mtx_def one_sq_mtx_def)
+  using assms by simp_all
+
+lemma mtx_invertible_inv[simp]: "mtx_invertible A \<Longrightarrow> mtx_invertible (A\<^sup>-\<^sup>1)"
+  using mtx_invertibleD mtx_invertibleI by blast
+
+lemma mtx_invertible_one[simp]: "mtx_invertible 1"
+  by (simp add: one_sq_mtx.rep_eq)
+
+lemma sq_mtx_inv_unique:
+  assumes "A * B = 1" and "B * A = 1"
+  shows "A\<^sup>-\<^sup>1 = B"
+  by (metis (no_types, lifting) assms mtx_invertibleD(2) 
+      mtx_invertibleI mult.assoc sq_mtx_one_idty(1))
+
+lemma sq_mtx_inv_idempotent[simp]: "mtx_invertible A \<Longrightarrow> A\<^sup>-\<^sup>1\<^sup>-\<^sup>1 = A"
+  using mtx_invertibleD sq_mtx_inv_unique by blast
+
+lemma sq_mtx_inv_mult:
+  assumes "mtx_invertible A" and "mtx_invertible B"
+  shows "(A * B)\<^sup>-\<^sup>1 = B\<^sup>-\<^sup>1 * A\<^sup>-\<^sup>1"
+  by (simp add: assms matrix_inv_matrix_mul sq_mtx_inv_def times_sq_mtx_def)
+
+lemma sq_mtx_inv_one[simp]: "1\<^sup>-\<^sup>1 = 1"
+  by (simp add: sq_mtx_inv_unique)
+
+definition similar_sq_mtx :: "('n::finite) sq_mtx \<Rightarrow> 'n sq_mtx \<Rightarrow> bool" (infixr "\<sim>" 25)
+  where "(A \<sim> B) \<longleftrightarrow> (\<exists> P. mtx_invertible P \<and> A = P\<^sup>-\<^sup>1 * B * P)"
+
+lemma similar_sq_mtx_matrix: "(A \<sim> B) = similar_matrix (to_vec A) (to_vec B)"
+  apply(unfold similar_matrix_def similar_sq_mtx_def, safe)
+   apply (metis sq_mtx_inv.rep_eq times_sq_mtx.rep_eq)
+  by (metis UNIV_I sq_mtx_inv.abs_eq times_sq_mtx.abs_eq to_mtx_inverse to_vec_inverse)
+
+lemma similar_sq_mtx_refl[simp]: "A \<sim> A"
+  by (unfold similar_sq_mtx_def, rule_tac x="1" in exI, simp)
+
+lemma similar_sq_mtx_simm: "A \<sim> B \<Longrightarrow> B \<sim> A"
+  apply(unfold similar_sq_mtx_def, clarsimp)
+  apply(rule_tac x="P\<^sup>-\<^sup>1" in exI, simp add: mult.assoc)
+  by (metis mtx_invertibleD(2) mult.assoc mult.left_neutral)
+
+lemma similar_sq_mtx_trans: "A \<sim> B \<Longrightarrow> B \<sim> C \<Longrightarrow> A \<sim> C"
+  unfolding similar_sq_mtx_matrix using similar_matrix_trans by blast
+
+lemma power_sq_mtx_diag: "(sq_mtx_diag f)^n = (\<d>\<i>\<a>\<g> i. f i^n)"
+  by (induct n, simp_all)
+
+lemma power_similiar_sq_mtx_diag_eq:
+  assumes "mtx_invertible P"
+      and "A = P\<^sup>-\<^sup>1 * (sq_mtx_diag f) * P"
+    shows "A^n = P\<^sup>-\<^sup>1 * (\<d>\<i>\<a>\<g> i. f i^n) * P"
+proof(induct n, simp_all add: assms)
+  fix n::nat
+  have "P\<^sup>-\<^sup>1 * sq_mtx_diag f * P * (P\<^sup>-\<^sup>1 * (\<d>\<i>\<a>\<g> i. f i ^ n) * P) = 
+  P\<^sup>-\<^sup>1 * sq_mtx_diag f * (\<d>\<i>\<a>\<g> i. f i ^ n) * P"
+    by (metis (no_types, lifting) assms(1) mtx_invertibleD(2) mult.assoc mult.right_neutral)
+  also have "... = P\<^sup>-\<^sup>1 * (\<d>\<i>\<a>\<g> i. f i * f i ^ n) * P"
+    by (simp add: mult.assoc) 
+  finally show "P\<^sup>-\<^sup>1 * sq_mtx_diag f * P * (P\<^sup>-\<^sup>1 * (\<d>\<i>\<a>\<g> i. f i ^ n) * P) = 
+  P\<^sup>-\<^sup>1 * (\<d>\<i>\<a>\<g> i. f i * f i ^ n) * P" .
+qed
+
+lemma power_similar_sq_mtx_diag:
+  assumes "A \<sim> (sq_mtx_diag f)"
+  shows "A^n \<sim> (\<d>\<i>\<a>\<g> i. f i^n)"
+  using assms power_similiar_sq_mtx_diag_eq 
+  unfolding similar_sq_mtx_def by blast
+
+
+subsection \<open> Banach space of square matrices \<close>
+
+lemma Cauchy_cols:
+  fixes X :: "nat \<Rightarrow> ('a::finite) sq_mtx" 
+  assumes "Cauchy X"
+  shows "Cauchy (\<lambda>n. \<c>\<o>\<l> i (X n))" 
+proof(unfold Cauchy_def dist_norm, clarsimp)
+  fix \<epsilon>::real assume "\<epsilon> > 0"
+  then obtain M where M_def:"\<forall>m\<ge>M. \<forall>n\<ge>M. \<parallel>X m - X n\<parallel> < \<epsilon>"
+    using \<open>Cauchy X\<close> unfolding Cauchy_def by(simp add: dist_sq_mtx_def) metis
+  {fix m n assume "m \<ge> M" and "n \<ge> M"
+    hence "\<epsilon> > \<parallel>X m - X n\<parallel>" 
+      using M_def by blast
+    moreover have "\<parallel>X m - X n\<parallel> \<ge> \<parallel>(X m - X n) *\<^sub>V \<e> i\<parallel>"
+      by(rule le_mtx_norm[of _ "X m - X n"], force)
+    moreover have "\<parallel>(X m - X n) *\<^sub>V \<e> i\<parallel> = \<parallel>X m *\<^sub>V \<e> i - X n *\<^sub>V \<e> i\<parallel>"
+      by (simp add: mtx_vec_mult_minus_rdistrib)
+    moreover have "... = \<parallel>\<c>\<o>\<l> i (X m) - \<c>\<o>\<l> i (X n)\<parallel>"
+      by (simp add: mtx_vec_mult_minus_rdistrib mtx_vec_mult_canon)
+    ultimately have "\<parallel>\<c>\<o>\<l> i (X m) - \<c>\<o>\<l> i (X n)\<parallel> < \<epsilon>" 
+      by linarith}
+  thus "\<exists>M. \<forall>m\<ge>M. \<forall>n\<ge>M. \<parallel>\<c>\<o>\<l> i (X m) - \<c>\<o>\<l> i (X n)\<parallel> < \<epsilon>" 
+    by blast
+qed
+
+lemma col_convergence:
+  assumes "\<forall>i. (\<lambda>n. \<c>\<o>\<l> i (X n)) \<longlonglongrightarrow> L $ i" 
+  shows "X \<longlonglongrightarrow> to_mtx (transpose L)"
+proof(unfold LIMSEQ_def dist_norm, clarsimp)
+  let ?L = "to_mtx (transpose L)"
+  let ?a = "CARD('a)" fix \<epsilon>::real assume "\<epsilon> > 0"
+  hence "\<epsilon> / ?a > 0" by simp
+  hence "\<forall>i. \<exists> N. \<forall>n\<ge>N. \<parallel>\<c>\<o>\<l> i (X n) - L $ i\<parallel> < \<epsilon>/?a"
+    using assms unfolding LIMSEQ_def dist_norm convergent_def by blast
+  then obtain N where "\<forall>i. \<forall>n\<ge>N. \<parallel>\<c>\<o>\<l> i (X n) - L $ i\<parallel> < \<epsilon>/?a"
+    using finite_nat_minimal_witness[of "\<lambda> i n. \<parallel>\<c>\<o>\<l> i (X n) - L $ i\<parallel> < \<epsilon>/?a"] by blast
+  also have "\<And>i n. (\<c>\<o>\<l> i (X n) - L $ i) = (\<c>\<o>\<l> i (X n - ?L))"
+    unfolding minus_sq_mtx_def by(transfer, simp add: transpose_def vec_eq_iff column_def)
+  ultimately have N_def:"\<forall>i. \<forall>n\<ge>N. \<parallel>\<c>\<o>\<l> i (X n - ?L)\<parallel> < \<epsilon>/?a" 
+    by auto
+  have "\<forall>n\<ge>N. \<parallel>X n - ?L\<parallel> < \<epsilon>"
+  proof(rule allI, rule impI)
+    fix n::nat assume "N \<le> n"
+    hence "\<forall> i. \<parallel>\<c>\<o>\<l> i (X n - ?L)\<parallel> < \<epsilon>/?a"
+      using N_def by blast
+    hence "(\<Sum>i\<in>UNIV. \<parallel>\<c>\<o>\<l> i (X n - ?L)\<parallel>) < (\<Sum>(i::'a)\<in>UNIV. \<epsilon>/?a)"
+      using sum_strict_mono[of _ "\<lambda>i. \<parallel>\<c>\<o>\<l> i (X n - ?L)\<parallel>"] by force
+    moreover have "\<parallel>X n - ?L\<parallel> \<le> (\<Sum>i\<in>UNIV. \<parallel>\<c>\<o>\<l> i (X n - ?L)\<parallel>)"
+      using sq_mtx_norm_le_sum_col by blast
+    moreover have "(\<Sum>(i::'a)\<in>UNIV. \<epsilon>/?a) = \<epsilon>" 
+      by force
+    ultimately show "\<parallel>X n - ?L\<parallel> < \<epsilon>" 
+      by linarith
+  qed
+  thus "\<exists>no. \<forall>n\<ge>no. \<parallel>X n - ?L\<parallel> < \<epsilon>" 
+    by blast
+qed
+
+instance sq_mtx :: (finite) banach
+proof(standard)
+  fix X :: "nat \<Rightarrow> 'a sq_mtx"
+  assume "Cauchy X"
+  hence "\<And>i. Cauchy (\<lambda>n. \<c>\<o>\<l> i (X n))"
+    using Cauchy_cols by blast
+  hence obs: "\<forall>i. \<exists>! L. (\<lambda>n. \<c>\<o>\<l> i (X n)) \<longlonglongrightarrow> L"
+    using Cauchy_convergent convergent_def LIMSEQ_unique by fastforce
+  define L where "L = (\<chi> i. lim (\<lambda>n. \<c>\<o>\<l> i (X n)))"
+  hence "\<forall>i. (\<lambda>n. \<c>\<o>\<l> i (X n)) \<longlonglongrightarrow> L $ i" 
+    using obs theI_unique[of "\<lambda>L. (\<lambda>n. \<c>\<o>\<l> _ (X n)) \<longlonglongrightarrow> L" "L $ _"] by (simp add: lim_def)
+  thus "convergent X"
+    using col_convergence unfolding convergent_def by blast
+qed
+
+lemma exp_similiar_sq_mtx_diag_eq:
+  assumes "mtx_invertible P"
+      and "A = P\<^sup>-\<^sup>1 * (\<d>\<i>\<a>\<g> i. f i) * P"
+    shows "exp A = P\<^sup>-\<^sup>1 * exp (\<d>\<i>\<a>\<g> i. f i) * P"
+proof(unfold exp_def power_similiar_sq_mtx_diag_eq[OF assms])
+  have "(\<Sum>n. P\<^sup>-\<^sup>1 * (\<d>\<i>\<a>\<g> i. f i ^ n) * P /\<^sub>R fact n) = 
+  (\<Sum>n. P\<^sup>-\<^sup>1 * ((\<d>\<i>\<a>\<g> i. f i ^ n) /\<^sub>R fact n) * P)"
+    by simp
+  also have "... = (\<Sum>n. P\<^sup>-\<^sup>1 * ((\<d>\<i>\<a>\<g> i. f i ^ n) /\<^sub>R fact n)) * P"
+    apply(subst suminf_multr[OF bounded_linear.summable[OF bounded_linear_mult_right]])
+    unfolding power_sq_mtx_diag[symmetric] by (simp_all add: summable_exp_generic)
+  also have "... = P\<^sup>-\<^sup>1 * (\<Sum>n. (\<d>\<i>\<a>\<g> i. f i ^ n) /\<^sub>R fact n) * P"
+    apply(subst suminf_mult[of _ "P\<^sup>-\<^sup>1"])
+    unfolding power_sq_mtx_diag[symmetric] 
+    by (simp_all add: summable_exp_generic)
+  finally show "(\<Sum>n. P\<^sup>-\<^sup>1 * (\<d>\<i>\<a>\<g> i. f i ^ n) * P /\<^sub>R fact n) = 
+  P\<^sup>-\<^sup>1 * (\<Sum>n. sq_mtx_diag f ^ n /\<^sub>R fact n) * P"
+    unfolding power_sq_mtx_diag by simp
+qed
+
+lemma exp_similiar_sq_mtx_diag:
+  assumes "A \<sim> sq_mtx_diag f"
+  shows "exp A \<sim> exp (sq_mtx_diag f)"
+  using assms exp_similiar_sq_mtx_diag_eq 
+  unfolding similar_sq_mtx_def by blast
+
+lemma suminf_sq_mtx_diag:
+  assumes "\<forall>i. (\<lambda>n. f n i) sums (suminf (\<lambda>n. f n i))"
+  shows "(\<Sum>n. (\<d>\<i>\<a>\<g> i. f n i)) = (\<d>\<i>\<a>\<g> i. \<Sum>n. f n i)"
+proof(rule suminfI, unfold sums_def LIMSEQ_iff, clarsimp simp: norm_sq_mtx_diag)
+  let ?g = "\<lambda>n i. \<bar>(\<Sum>n<n. f n i) - (\<Sum>n. f n i)\<bar>"
+  fix r::real assume "r > 0"
+  have "\<forall>i. \<exists>no. \<forall>n\<ge>no. ?g n i < r"
+    using assms \<open>r > 0\<close> unfolding sums_def LIMSEQ_iff by clarsimp 
+  then obtain N where key: "\<forall>i. \<forall>n\<ge>N. ?g n i < r"
+    using finite_nat_minimal_witness[of "\<lambda>i n. ?g n i < r"] by blast
+  {fix n::nat
+    assume "n \<ge> N"
+    obtain i where i_def: "Max {x. \<exists>i. x = ?g n i} = ?g n i"
+      using cMax_finite_ex[of "{x. \<exists>i. x = ?g n i}"] by auto
+    hence "?g n i < r"
+      using key \<open>n \<ge> N\<close> by blast
+    hence "Max {x. \<exists>i. x = ?g n i} < r"
+      unfolding i_def[symmetric] .}
+  thus "\<exists>N. \<forall>n\<ge>N. Max {x. \<exists>i. x = ?g n i} < r"
+    by blast
+qed
+
+lemma exp_sq_mtx_diag: "exp (sq_mtx_diag f) = (\<d>\<i>\<a>\<g> i. exp (f i))"
+  apply(unfold exp_def, simp add: power_sq_mtx_diag scaleR_sq_mtx_diag)
+  apply(rule suminf_sq_mtx_diag)
+  using exp_converges[of "f _"] 
+  unfolding sums_def LIMSEQ_iff exp_def by force
+
+lemma exp_scaleR_diagonal1:
+  assumes "mtx_invertible P" and "A = P\<^sup>-\<^sup>1 * (\<d>\<i>\<a>\<g> i. f i) * P"
+    shows "exp (t *\<^sub>R A) = P\<^sup>-\<^sup>1 * (\<d>\<i>\<a>\<g> i. exp (t * f i)) * P"
+proof-
+  have "exp (t *\<^sub>R A) = exp (P\<^sup>-\<^sup>1 * (t *\<^sub>R sq_mtx_diag f) * P)"
+    using assms by simp
+  also have "... = P\<^sup>-\<^sup>1 * (\<d>\<i>\<a>\<g> i. exp (t * f i)) * P"
+    by (metis assms(1) exp_similiar_sq_mtx_diag_eq exp_sq_mtx_diag scaleR_sq_mtx_diag)
+  finally show "exp (t *\<^sub>R A) = P\<^sup>-\<^sup>1 * (\<d>\<i>\<a>\<g> i. exp (t * f i)) * P" .
+qed
+
+lemma exp_scaleR_diagonal2:
+  assumes "mtx_invertible P" and "A = P * (\<d>\<i>\<a>\<g> i. f i) * P\<^sup>-\<^sup>1"
+    shows "exp (t *\<^sub>R A) = P * (\<d>\<i>\<a>\<g> i. exp (t * f i)) * P\<^sup>-\<^sup>1"
+  apply(subst sq_mtx_inv_idempotent[OF assms(1), symmetric])
+  apply(rule exp_scaleR_diagonal1)
+  by (simp_all add: assms)
+
+
+subsection \<open> Examples \<close>
+
+definition "mtx A = to_mtx (vector (map vector A))"
+
+lemma vector_nth_eq: "(vector A) $ i = foldr (\<lambda>x f n. (f (n + 1))(n := x)) A (\<lambda>n x. 0) 1 i"
+  unfolding vector_def by simp
+
+lemma mtx_ith_eq[simp]: "mtx A $$ i $ j = foldr (\<lambda>x f n. (f (n + 1))(n := x))
+  (map (\<lambda>l. vec_lambda (foldr (\<lambda>x f n. (f (n + 1))(n := x)) l (\<lambda>n x. 0) 1)) A) (\<lambda>n x. 0) 1 i $ j"
+  unfolding mtx_def vector_def by (simp add: vector_nth_eq)
+
+subsubsection \<open> 2x2 matrices \<close>
+
+lemma mtx2_eq_iff: "(mtx 
+  ([a1, b1] # 
+   [c1, d1] # []) :: 2 sq_mtx) = mtx 
+  ([a2, b2] # 
+   [c2, d2] # []) \<longleftrightarrow> a1 = a2 \<and> b1 = b2 \<and> c1 = c2 \<and> d1 = d2"
+  apply(simp add: sq_mtx_eq_iff, safe)
+  using exhaust_2 by force+
+
+lemma mtx2_to_mtx: "mtx 
+  ([a, b] # 
+   [c, d] # []) = 
+  to_mtx (\<chi> i j::2. if i=1 \<and> j=1 then a 
+  else (if i=1 \<and> j=2 then b 
+  else (if i=2 \<and> j=1 then c 
+  else d)))"
+  apply(subst sq_mtx_eq_iff)
+  using exhaust_2 by force
+
+abbreviation diag2 :: "real \<Rightarrow> real \<Rightarrow> 2 sq_mtx" 
+  where "diag2 \<iota>\<^sub>1 \<iota>\<^sub>2 \<equiv> mtx 
+   ([\<iota>\<^sub>1, 0] # 
+    [0, \<iota>\<^sub>2] # [])"
+
+lemma diag2_eq: "diag2 (\<iota> 1) (\<iota> 2) = (\<d>\<i>\<a>\<g> i. \<iota> i)"
+  apply(simp add: sq_mtx_eq_iff)
+  using exhaust_2 by (force simp: axis_def)
+
+lemma one_mtx2: "(1::2 sq_mtx) = diag2 1 1"
+  apply(subst sq_mtx_eq_iff)
+  using exhaust_2 by force
+
+lemma zero_mtx2: "(0::2 sq_mtx) = diag2 0 0"
+  by (simp add: sq_mtx_eq_iff)
+
+lemma scaleR_mtx2: "k *\<^sub>R mtx 
+  ([a, b] # 
+   [c, d] # []) = mtx 
+  ([k*a, k*b] # 
+   [k*c, k*d] # [])"
+  by (simp add: sq_mtx_eq_iff)
+
+lemma uminus_mtx2: "-mtx 
+  ([a, b] # 
+   [c, d] # []) = (mtx 
+  ([-a, -b] # 
+   [-c, -d] # [])::2 sq_mtx)"
+  by (simp add: sq_mtx_uminus_eq sq_mtx_eq_iff)
+
+lemma plus_mtx2: "mtx 
+  ([a1, b1] # 
+   [c1, d1] # []) + mtx 
+  ([a2, b2] # 
+   [c2, d2] # []) = ((mtx 
+  ([a1+a2, b1+b2] # 
+   [c1+c2, d1+d2] # []))::2 sq_mtx)"
+  by (simp add: sq_mtx_eq_iff)
+
+lemma minus_mtx2: "mtx 
+  ([a1, b1] # 
+   [c1, d1] # []) - mtx 
+  ([a2, b2] # 
+   [c2, d2] # []) = ((mtx 
+  ([a1-a2, b1-b2] # 
+   [c1-c2, d1-d2] # []))::2 sq_mtx)"
+  by (simp add: sq_mtx_eq_iff)
+
+lemma times_mtx2: "mtx 
+  ([a1, b1] # 
+   [c1, d1] # []) * mtx 
+  ([a2, b2] # 
+   [c2, d2] # []) = ((mtx 
+  ([a1*a2+b1*c2, a1*b2+b1*d2] # 
+   [c1*a2+d1*c2, c1*b2+d1*d2] # []))::2 sq_mtx)"
+  unfolding sq_mtx_times_eq UNIV_2
+  by (simp add: sq_mtx_eq_iff)
+
+subsubsection \<open> 3x3 matrices \<close>
+
+lemma mtx3_to_mtx: "mtx 
+  ([a\<^sub>1\<^sub>1, a\<^sub>1\<^sub>2, a\<^sub>1\<^sub>3] # 
+   [a\<^sub>2\<^sub>1, a\<^sub>2\<^sub>2, a\<^sub>2\<^sub>3] # 
+   [a\<^sub>3\<^sub>1, a\<^sub>3\<^sub>2, a\<^sub>3\<^sub>3] # []) = 
+  to_mtx (\<chi> i j::3. if i=1 \<and> j=1 then a\<^sub>1\<^sub>1
+  else (if i=1 \<and> j=2 then a\<^sub>1\<^sub>2 
+  else (if i=1 \<and> j=3 then a\<^sub>1\<^sub>3 
+  else (if i=2 \<and> j=1 then a\<^sub>2\<^sub>1
+  else (if i=2 \<and> j=2 then a\<^sub>2\<^sub>2 
+  else (if i=2 \<and> j=3 then a\<^sub>2\<^sub>3 
+  else (if i=3 \<and> j=1 then a\<^sub>3\<^sub>1 
+  else (if i=3 \<and> j=2 then a\<^sub>3\<^sub>2 
+  else a\<^sub>3\<^sub>3))))))))"
+  apply(simp add: sq_mtx_eq_iff)
+  using exhaust_3 by force
+
+abbreviation diag3 :: "real \<Rightarrow> real \<Rightarrow> real \<Rightarrow> 3 sq_mtx" 
+  where "diag3 \<iota>\<^sub>1 \<iota>\<^sub>2 \<iota>\<^sub>3 \<equiv> mtx 
+  ([\<iota>\<^sub>1, 0, 0] # 
+   [0, \<iota>\<^sub>2, 0] # 
+   [0, 0, \<iota>\<^sub>3] # [])"
+
+lemma diag3_eq: "diag3 (\<iota> 1) (\<iota> 2) (\<iota> 3) = (\<d>\<i>\<a>\<g> i. \<iota> i)"
+  apply(simp add: sq_mtx_eq_iff)
+  using exhaust_3 by (force simp: axis_def)
+
+lemma one_mtx3: "(1::3 sq_mtx) = diag3 1 1 1"
+  apply(subst sq_mtx_eq_iff)
+  using exhaust_3 by force
+
+lemma zero_mtx3: "(0::3 sq_mtx) = diag3 0 0 0"
+  by (simp add: sq_mtx_eq_iff)
+
+lemma scaleR_mtx3: "k *\<^sub>R mtx 
+  ([a\<^sub>1\<^sub>1, a\<^sub>1\<^sub>2, a\<^sub>1\<^sub>3] # 
+   [a\<^sub>2\<^sub>1, a\<^sub>2\<^sub>2, a\<^sub>2\<^sub>3] # 
+   [a\<^sub>3\<^sub>1, a\<^sub>3\<^sub>2, a\<^sub>3\<^sub>3] # []) = mtx 
+  ([k*a\<^sub>1\<^sub>1, k*a\<^sub>1\<^sub>2, k*a\<^sub>1\<^sub>3] # 
+   [k*a\<^sub>2\<^sub>1, k*a\<^sub>2\<^sub>2, k*a\<^sub>2\<^sub>3] # 
+   [k*a\<^sub>3\<^sub>1, k*a\<^sub>3\<^sub>2, k*a\<^sub>3\<^sub>3] # [])"
+  by (simp add: sq_mtx_eq_iff)
+
+lemma plus_mtx3: "mtx 
+  ([a\<^sub>1\<^sub>1, a\<^sub>1\<^sub>2, a\<^sub>1\<^sub>3] # 
+   [a\<^sub>2\<^sub>1, a\<^sub>2\<^sub>2, a\<^sub>2\<^sub>3] # 
+   [a\<^sub>3\<^sub>1, a\<^sub>3\<^sub>2, a\<^sub>3\<^sub>3] # []) + mtx 
+  ([b\<^sub>1\<^sub>1, b\<^sub>1\<^sub>2, b\<^sub>1\<^sub>3] # 
+   [b\<^sub>2\<^sub>1, b\<^sub>2\<^sub>2, b\<^sub>2\<^sub>3] # 
+   [b\<^sub>3\<^sub>1, b\<^sub>3\<^sub>2, b\<^sub>3\<^sub>3] # []) = (mtx 
+  ([a\<^sub>1\<^sub>1+b\<^sub>1\<^sub>1, a\<^sub>1\<^sub>2+b\<^sub>1\<^sub>2, a\<^sub>1\<^sub>3+b\<^sub>1\<^sub>3] # 
+   [a\<^sub>2\<^sub>1+b\<^sub>2\<^sub>1, a\<^sub>2\<^sub>2+b\<^sub>2\<^sub>2, a\<^sub>2\<^sub>3+b\<^sub>2\<^sub>3] # 
+   [a\<^sub>3\<^sub>1+b\<^sub>3\<^sub>1, a\<^sub>3\<^sub>2+b\<^sub>3\<^sub>2, a\<^sub>3\<^sub>3+b\<^sub>3\<^sub>3] # [])::3 sq_mtx)"
+  by (subst sq_mtx_eq_iff) simp
+
+lemma minus_mtx3: "mtx 
+  ([a\<^sub>1\<^sub>1, a\<^sub>1\<^sub>2, a\<^sub>1\<^sub>3] # 
+   [a\<^sub>2\<^sub>1, a\<^sub>2\<^sub>2, a\<^sub>2\<^sub>3] # 
+   [a\<^sub>3\<^sub>1, a\<^sub>3\<^sub>2, a\<^sub>3\<^sub>3] # []) - mtx 
+  ([b\<^sub>1\<^sub>1, b\<^sub>1\<^sub>2, b\<^sub>1\<^sub>3] # 
+   [b\<^sub>2\<^sub>1, b\<^sub>2\<^sub>2, b\<^sub>2\<^sub>3] # 
+   [b\<^sub>3\<^sub>1, b\<^sub>3\<^sub>2, b\<^sub>3\<^sub>3] # []) = (mtx 
+  ([a\<^sub>1\<^sub>1-b\<^sub>1\<^sub>1, a\<^sub>1\<^sub>2-b\<^sub>1\<^sub>2, a\<^sub>1\<^sub>3-b\<^sub>1\<^sub>3] # 
+   [a\<^sub>2\<^sub>1-b\<^sub>2\<^sub>1, a\<^sub>2\<^sub>2-b\<^sub>2\<^sub>2, a\<^sub>2\<^sub>3-b\<^sub>2\<^sub>3] # 
+   [a\<^sub>3\<^sub>1-b\<^sub>3\<^sub>1, a\<^sub>3\<^sub>2-b\<^sub>3\<^sub>2, a\<^sub>3\<^sub>3-b\<^sub>3\<^sub>3] # [])::3 sq_mtx)"
+  by (simp add: sq_mtx_eq_iff)
+
+lemma times_mtx3: "mtx 
+  ([a\<^sub>1\<^sub>1, a\<^sub>1\<^sub>2, a\<^sub>1\<^sub>3] # 
+   [a\<^sub>2\<^sub>1, a\<^sub>2\<^sub>2, a\<^sub>2\<^sub>3] # 
+   [a\<^sub>3\<^sub>1, a\<^sub>3\<^sub>2, a\<^sub>3\<^sub>3] # []) * mtx 
+  ([b\<^sub>1\<^sub>1, b\<^sub>1\<^sub>2, b\<^sub>1\<^sub>3] # 
+   [b\<^sub>2\<^sub>1, b\<^sub>2\<^sub>2, b\<^sub>2\<^sub>3] # 
+   [b\<^sub>3\<^sub>1, b\<^sub>3\<^sub>2, b\<^sub>3\<^sub>3] # []) = (mtx 
+  ([a\<^sub>1\<^sub>1*b\<^sub>1\<^sub>1+a\<^sub>1\<^sub>2*b\<^sub>2\<^sub>1+a\<^sub>1\<^sub>3*b\<^sub>3\<^sub>1, a\<^sub>1\<^sub>1*b\<^sub>1\<^sub>2+a\<^sub>1\<^sub>2*b\<^sub>2\<^sub>2+a\<^sub>1\<^sub>3*b\<^sub>3\<^sub>2, a\<^sub>1\<^sub>1*b\<^sub>1\<^sub>3+a\<^sub>1\<^sub>2*b\<^sub>2\<^sub>3+a\<^sub>1\<^sub>3*b\<^sub>3\<^sub>3] # 
+   [a\<^sub>2\<^sub>1*b\<^sub>1\<^sub>1+a\<^sub>2\<^sub>2*b\<^sub>2\<^sub>1+a\<^sub>2\<^sub>3*b\<^sub>3\<^sub>1, a\<^sub>2\<^sub>1*b\<^sub>1\<^sub>2+a\<^sub>2\<^sub>2*b\<^sub>2\<^sub>2+a\<^sub>2\<^sub>3*b\<^sub>3\<^sub>2, a\<^sub>2\<^sub>1*b\<^sub>1\<^sub>3+a\<^sub>2\<^sub>2*b\<^sub>2\<^sub>3+a\<^sub>2\<^sub>3*b\<^sub>3\<^sub>3] # 
+   [a\<^sub>3\<^sub>1*b\<^sub>1\<^sub>1+a\<^sub>3\<^sub>2*b\<^sub>2\<^sub>1+a\<^sub>3\<^sub>3*b\<^sub>3\<^sub>1, a\<^sub>3\<^sub>1*b\<^sub>1\<^sub>2+a\<^sub>3\<^sub>2*b\<^sub>2\<^sub>2+a\<^sub>3\<^sub>3*b\<^sub>3\<^sub>2, a\<^sub>3\<^sub>1*b\<^sub>1\<^sub>3+a\<^sub>3\<^sub>2*b\<^sub>2\<^sub>3+a\<^sub>3\<^sub>3*b\<^sub>3\<^sub>3] # [])::3 sq_mtx)"
+  unfolding sq_mtx_times_eq
+  unfolding UNIV_3 by (simp add: sq_mtx_eq_iff)
+
+end
\ No newline at end of file
diff --git a/thys/Matrices_for_ODEs/document/root.bib b/thys/Matrices_for_ODEs/document/root.bib
new file mode 100644
--- /dev/null
+++ b/thys/Matrices_for_ODEs/document/root.bib
@@ -0,0 +1,70 @@
+%% This BibTeX bibliography file was created using BibDesk.
+%% http://bibdesk.sourceforge.net/
+
+
+%% Created for Jonathan Julian Huerta y Munive at 2019-06-10 16:26:22 +0100 
+
+
+%% Saved with string encoding Unicode (UTF-8)
+
+@article{ArmstrongGS16,
+	Author = {Alasdair Armstrong and Victor B. F. Gomes and Georg Struth},
+	Journal = {Formal Aspects of Computing},
+	Number = 2,
+	Pages = {265--293},
+	Title = {Building program construction and verification tools from algebraic principles},
+	Volume = 28,
+	Year = 2016}
+
+@Unpublished{FosterMS19,
+  author    = {Simon Foster and
+               Jonathan Juli{\'{a}}n Huerta y Munive and
+               Georg Struth},
+  title     = {Differential {H}oare Logics and Refinement Calculi for Hybrid Systems
+               with {I}sabelle/{HOL}},
+  journal   = {CoRR},
+  volume    = {abs/1910.13554},
+  note = {\href{http://arxiv.org/abs/1909.05618}{arXiv:1909.05618}[cs.LO]},
+  year      = {2019}
+}
+
+@Unpublished{MuniveS19,
+  author = 	 {Huerta y Munive, Jonathan Juli\'an and Struth, G.},
+  title = 	 {Predicate Transformer Semantics for Hybrid Systems: Verification Components for {I}sabelle/{HOL}},
+    note =         {\href{https://arxiv.org/abs/arXiv:1909.05618}{arXiv:1909.05618} [cs.LO]},
+      year =         2019,
+}
+
+@article{ImmlerH12a,
+	Author = {Fabian Immler and Johannes H{\"{o}}lzl},
+	Journal = {Archive of Formal Proofs},
+	Title = {Ordinary Differential Equations},
+	Url = {https://www.isa-afp.org/entries/Ordinary_Differential_Equations.shtml},
+	Year = {2012},
+	Bdsk-Url-1 = {https://www.isa-afp.org/entries/Ordinary_Differential_Equations.shtml}}
+
+@book{Platzer10,
+	Author = {Andr\'e Platzer},
+	Publisher = {Springer},
+	Title = {Logical Analysis of Hybrid Systems},
+	Year = 2010}
+
+@article{afp:hybrid,
+	Author = {Huerta y Munive, Jonathan Juli\'an},
+	Journal = {Archive of Formal Proofs},
+	Title = {Verification Components for Hybrid Systems},
+	Year = 2019}
+
+@article{afp:transem,
+	Author = {Georg Struth},
+	Journal = {Archive of Formal Proofs},
+	Title = {Transformer Semantics},
+	Year = 2018}
+
+@article{afp:vericomp,
+	Author = {Victor B. F. Gomes and Georg Struth},
+	Journal = {Archive of Formal Proofs},
+	Title = {Program Construction and Verification Components Based on {K}leene Algebra},
+	Year = 2016}
+
+
diff --git a/thys/Matrices_for_ODEs/document/root.tex b/thys/Matrices_for_ODEs/document/root.tex
new file mode 100644
--- /dev/null
+++ b/thys/Matrices_for_ODEs/document/root.tex
@@ -0,0 +1,76 @@
+\documentclass[11pt,a4paper]{article}
+\usepackage{isabelle,isabellesym}
+
+% further packages required for unusual symbols (see also
+% isabellesym.sty), use only when needed
+
+\usepackage{amssymb}
+  %for \<leadsto>, \<box>, \<diamond>, \<sqsupset>, \<mho>, \<Join>,
+  %\<lhd>, \<lesssim>, \<greatersim>, \<lessapprox>, \<greaterapprox>,
+  %\<triangleq>, \<yen>, \<lozenge>
+
+%\usepackage{eurosym}
+  %for \<euro>
+
+%\usepackage[only,bigsqcap]{stmaryrd}
+  %for \<Sqinter>
+
+%\usepackage{eufrak}
+  %for \<AA> ... \<ZZ>, \<aa> ... \<zz> (also included in amssymb)
+
+%\usepackage{textcomp}
+  %for \<onequarter>, \<onehalf>, \<threequarters>, \<degree>, \<cent>,
+  %\<currency>
+
+% this should be the last package used
+\usepackage{pdfsetup}
+
+% urls in roman style, theory text in math-similar italics
+\urlstyle{rm}
+\isabellestyle{it}
+
+% for uniform font size
+%\renewcommand{\isastyle}{\isastyleminor}
+
+\renewcommand{\isasymlonglonglongrightarrow}{$\longrightarrow$}
+
+
+\begin{document}
+
+\title{Matrices for ODEs}
+\author{Jonathan Juli\'an Huerta y Munive}
+\maketitle
+
+\begin{abstract}
+  Our theories formalise various matrix properties that serve to establish 
+  existence, uniqueness and characterisation of the solution to affine 
+  systems of ordinary differential equations (ODEs). In particular, we
+  formalise the operator and maximum norm of matrices. Then we use 
+  them to prove that square matrices form a Banach space, and in this 
+  setting, we show an instance of Picard-Lindel\"of's theorem for affine
+  systems of ODEs. Finally, we apply this formalisation by verifying three 
+  simple hybrid programs.
+\end{abstract}
+
+\tableofcontents
+
+% sane default for proof documents
+\parindent 0pt\parskip 0.5ex
+
+\section{Introductory Remarks}
+
+Affine systems of ordinary differential equations (ODEs) are those whose associated vector fields are linear transformations. That is, if there is a matrix-valued function $A:\mathbb{R}\to M_{n\times n}(\mathbb{R})$ and vector function $B:\mathbb{R}\to\mathbb{R}^n$ such that the system of ODEs $x'\, t=f\, (t,x\, t)$ can be rewritten as $x'\, t=A\cdot (x\, t)+B\, t$, then the system is affine. Similarly, the associated linear system of ODEs is $x'\, t=A\cdot (x\, t)$ for matrix-vector multiplication $\cdot$. Our theories formalise affine (hence linear) systems of ordinary differential equations. For this purpose, we extend the ODE libraries of~\cite{ImmlerH12a} and linear algebra in HOL-Analysis. We add to them various results about invertibility of matrices, their diagonalisation, their operator and maximum norms, and properties relating them with vectors. We also define a new type of square matrices and prove that this is a Banach space. Then we obtain results about derivatives of matrix-vector multiplication and use them to prove Picard-Lindel\"of's theorem as formalised in~\cite{afp:hybrid}. The Banach space instance allows us to characterise the general solution to affine systems of ODEs in terms of the matrix-exponential. Finally, we use the components of~\cite{afp:hybrid} to do three simple verification examples in the style of differential dynamic logic~\cite{Platzer10} as showcased in~\cite{ArmstrongGS16,FosterMS19,MuniveS19}. A paper with a detailed overview of the various contributions that this formalisation adds to the verification components will be available soon in ArXiv.
+
+% generated text of all theories
+\input{session}
+
+% optional bibliography
+\bibliographystyle{abbrv}
+\bibliography{root}
+
+\end{document}
+
+%%% Local Variables:
+%%% mode: latex
+%%% TeX-master: t
+%%% End:
diff --git a/thys/ROOTS b/thys/ROOTS
--- a/thys/ROOTS
+++ b/thys/ROOTS
@@ -1,537 +1,538 @@
 ADS_Functor
 AODV
 Attack_Trees
 Auto2_HOL
 Auto2_Imperative_HOL
 AVL-Trees
 AWN
 Abortable_Linearizable_Modules
 Abs_Int_ITP2012
 Abstract-Hoare-Logics
 Abstract-Rewriting
 Abstract_Completeness
 Abstract_Soundness
 Adaptive_State_Counting
 Affine_Arithmetic
 Aggregation_Algebras
 Akra_Bazzi
 Algebraic_Numbers
 Algebraic_VCs
 Allen_Calculus
 Amortized_Complexity
 AnselmGod
 Applicative_Lifting
 Approximation_Algorithms
 Architectural_Design_Patterns
 Aristotles_Assertoric_Syllogistic
 Arith_Prog_Rel_Primes
 ArrowImpossibilityGS
 AutoFocus-Stream
 Automatic_Refinement
 AxiomaticCategoryTheory
 BDD
 BNF_Operations
 Banach_Steinhaus
 Bell_Numbers_Spivey
 Berlekamp_Zassenhaus
 Bernoulli
 Bertrands_Postulate
 Bicategory
 BinarySearchTree
 Binding_Syntax_Theory
 Binomial-Heaps
 Binomial-Queues
 BNF_CC
 Bondy
 Boolean_Expression_Checkers
 Bounded_Deducibility_Security
 Buchi_Complementation
 Budan_Fourier
 Buffons_Needle
 Buildings
 BytecodeLogicJmlTypes
 C2KA_DistributedSystems
 CAVA_Automata
 CAVA_LTL_Modelchecker
 CCS
 CISC-Kernel
 CRDT
 CYK
 CakeML
 CakeML_Codegen
 Call_Arity
 Card_Equiv_Relations
 Card_Multisets
 Card_Number_Partitions
 Card_Partitions
 Cartan_FP
 Case_Labeling
 Catalan_Numbers
 Category
 Category2
 Category3
 Cauchy
 Cayley_Hamilton
 Certification_Monads
 Chord_Segments
 Circus
 Clean
 ClockSynchInst
 Closest_Pair_Points
 CofGroups
 Coinductive
 Coinductive_Languages
 Collections
 Comparison_Sort_Lower_Bound
 Compiling-Exceptions-Correctly
 Completeness
 Complete_Non_Orders
 Complex_Geometry
 Complx
 ComponentDependencies
 ConcurrentGC
 ConcurrentIMP
 Concurrent_Ref_Alg
 Concurrent_Revisions
 Consensus_Refined
 Constructive_Cryptography
 Constructor_Funs
 Containers
 CoreC++
 Core_DOM
 Count_Complex_Roots
 CryptHOL
 CryptoBasedCompositionalProperties
 DFS_Framework
 DPT-SAT-Solver
 DataRefinementIBP
 Datatype_Order_Generator
 Decl_Sem_Fun_PL
 Decreasing-Diagrams
 Decreasing-Diagrams-II
 Deep_Learning
 Density_Compiler
 Dependent_SIFUM_Refinement
 Dependent_SIFUM_Type_Systems
 Depth-First-Search
 Derangements
 Deriving
 Descartes_Sign_Rule
 Dict_Construction
 Differential_Dynamic_Logic
 Differential_Game_Logic
 Dijkstra_Shortest_Path
 Diophantine_Eqns_Lin_Hom
 Dirichlet_L
 Dirichlet_Series
 Discrete_Summation
 DiscretePricing
 DiskPaxos
 DynamicArchitectures
 Dynamic_Tables
 E_Transcendental
 Echelon_Form
 EdmondsKarp_Maxflow
 Efficient-Mergesort
 Elliptic_Curves_Group_Law
 Encodability_Process_Calculi
 Epistemic_Logic
 Ergodic_Theory
 Error_Function
 Euler_MacLaurin
 Euler_Partition
 Example-Submission
 Factored_Transition_System_Bounding
 Farkas
 FFT
 FLP
 FOL-Fitting
 FOL_Harrison
 FOL_Seq_Calc1
 Falling_Factorial_Sum
 FeatherweightJava
 Featherweight_OCL
 Fermat3_4
 FileRefinement
 FinFun
 Finger-Trees
 Finite_Automata_HF
 First_Order_Terms
 First_Welfare_Theorem
 Fishburn_Impossibility
 Fisher_Yates
 Flow_Networks
 Floyd_Warshall
 Flyspeck-Tame
 FocusStreamsCaseStudies
 Forcing
 Formal_SSA
 Formula_Derivatives
 Fourier
 Free-Boolean-Algebra
 Free-Groups
 FunWithFunctions
 FunWithTilings
 Functional-Automata
 Functional_Ordered_Resolution_Prover
 Furstenberg_Topology
 GPU_Kernel_PL
 Gabow_SCC
 Game_Based_Crypto
 Gauss-Jordan-Elim-Fun
 Gauss_Jordan
 Gauss_Sums
 Gaussian_Integers
 GenClock
 General-Triangle
 Generalized_Counting_Sort
 Generic_Deriving
 Generic_Join
 GewirthPGCProof
 Girth_Chromatic
 GoedelGod
 Goodstein_Lambda
 GraphMarkingIBP
 Graph_Saturation
 Graph_Theory
 Green
 Groebner_Bases
 Groebner_Macaulay
 Gromov_Hyperbolicity
 Group-Ring-Module
 HOL-CSP
 HOLCF-Prelude
 HRB-Slicing
 Heard_Of
 Hello_World
 HereditarilyFinite
 Hermite
 Hidden_Markov_Models
 Higher_Order_Terms
 Hoare_Time
 HotelKeyCards
 Huffman
 Hybrid_Logic
 Hybrid_Multi_Lane_Spatial_Logic
 Hybrid_Systems_VCs
 HyperCTL
 IEEE_Floating_Point
 IMAP-CRDT
 IMO2019
 IMP2
 IMP2_Binary_Heap
 IP_Addresses
 Imperative_Insertion_Sort
 Impossible_Geometry
 Incompleteness
 Incredible_Proof_Machine
 Inductive_Confidentiality
 InfPathElimination
 InformationFlowSlicing
 InformationFlowSlicing_Inter
 Integration
 Interval_Arithmetic_Word32
 Iptables_Semantics
 Irrationality_J_Hancl
 Isabelle_C
 Isabelle_Meta_Model
 Jacobson_Basic_Algebra
 Jinja
 JinjaThreads
 JiveDataStoreModel
 Jordan_Hoelder
 Jordan_Normal_Form
 KAD
 KAT_and_DRA
 KBPs
 KD_Tree
 Key_Agreement_Strong_Adversaries
 Kleene_Algebra
 Knot_Theory
 Knuth_Morris_Pratt
 Koenigsberg_Friendship
 Kruskal
 Kuratowski_Closure_Complement
 LLL_Basis_Reduction
 LLL_Factorization
 LOFT
 LTL
 LTL_to_DRA
 LTL_to_GBA
 LTL_Master_Theorem
 LTL_Normal_Form
 Lam-ml-Normalization
 LambdaAuth
 LambdaMu
 Lambda_Free_KBOs
 Lambda_Free_RPOs
 Lambert_W
 Landau_Symbols
 Laplace_Transform
 Latin_Square
 LatticeProperties
 Lambda_Free_EPO
 Launchbury
 Lazy-Lists-II
 Lazy_Case
 Lehmer
 Lifting_Definition_Option
 LightweightJava
 LinearQuantifierElim
 Linear_Inequalities
 Linear_Programming
 Linear_Recurrences
 Liouville_Numbers
 List-Index
 List-Infinite
 List_Interleaving
 List_Inversions
 List_Update
 LocalLexing
 Localization_Ring
 Locally-Nameless-Sigma
 Lowe_Ontological_Argument
 Lower_Semicontinuous
 Lp
 Lucas_Theorem
 MFMC_Countable
 MSO_Regex_Equivalence
 Markov_Models
 Marriage
 Mason_Stothers
+Matrices_for_ODEs
 Matrix
 Matrix_Tensor
 Matroids
 Max-Card-Matching
 Median_Of_Medians_Selection
 Menger
 Mersenne_Primes
 MFODL_Monitor_Optimized
 MFOTL_Monitor
 MiniML
 Minimal_SSA
 Minkowskis_Theorem
 Minsky_Machines
 Modal_Logics_for_NTS
 Modular_Assembly_Kit_Security
 Monad_Memo_DP
 Monad_Normalisation
 MonoBoolTranAlgebra
 MonoidalCategory
 Monomorphic_Monad
 MuchAdoAboutTwo
 Multirelations
 Multi_Party_Computation
 Myhill-Nerode
 Name_Carrying_Type_Inference
 Nat-Interval-Logic
 Native_Word
 Nested_Multisets_Ordinals
 Network_Security_Policy_Verification
 Neumann_Morgenstern_Utility
 No_FTL_observers
 Nominal2
 Noninterference_CSP
 Noninterference_Concurrent_Composition
 Noninterference_Generic_Unwinding
 Noninterference_Inductive_Unwinding
 Noninterference_Ipurge_Unwinding
 Noninterference_Sequential_Composition
 NormByEval
 Nullstellensatz
 Octonions
 Open_Induction
 OpSets
 Optics
 Optimal_BST
 Orbit_Stabiliser
 Order_Lattice_Props
 Ordered_Resolution_Prover
 Ordinal
 Ordinals_and_Cardinals
 Ordinary_Differential_Equations
 PCF
 PLM
 Pell
 POPLmark-deBruijn
 PSemigroupsConvolution
 Pairing_Heap
 Paraconsistency
 Parity_Game
 Partial_Function_MR
 Partial_Order_Reduction
 Password_Authentication_Protocol
 Perfect-Number-Thm
 Perron_Frobenius
 Pi_Calculus
 Pi_Transcendental
 Planarity_Certificates
 Polynomial_Factorization
 Polynomial_Interpolation
 Polynomials
 Poincare_Bendixson
 Poincare_Disc
 Pop_Refinement
 Posix-Lexing
 Possibilistic_Noninterference
 Power_Sum_Polynomials
 Pratt_Certificate
 Presburger-Automata
 Prim_Dijkstra_Simple
 Prime_Distribution_Elementary
 Prime_Harmonic_Series
 Prime_Number_Theorem
 Priority_Queue_Braun
 Priority_Search_Trees
 Probabilistic_Noninterference
 Probabilistic_Prime_Tests
 Probabilistic_System_Zoo
 Probabilistic_Timed_Automata
 Probabilistic_While
 Projective_Geometry
 Program-Conflict-Analysis
 Promela
 Proof_Strategy_Language
 PropResPI
 Propositional_Proof_Systems
 Prpu_Maxflow
 PseudoHoops
 Psi_Calculi
 Ptolemys_Theorem
 QHLProver
 QR_Decomposition
 Quantales
 Quaternions
 Quick_Sort_Cost
 RIPEMD-160-SPARK
 ROBDD
 RSAPSS
 Ramsey-Infinite
 Random_BSTs
 Randomised_BSTs
 Random_Graph_Subgraph_Threshold
 Randomised_Social_Choice
 Rank_Nullity_Theorem
 Real_Impl
 Recursion-Theory-I
 Refine_Imperative_HOL
 Refine_Monadic
 RefinementReactive
 Regex_Equivalence
 Regular-Sets
 Regular_Algebras
 Relation_Algebra
 Relational-Incorrectness-Logic
 Rep_Fin_Groups
 Residuated_Lattices
 Resolution_FOL
 Rewriting_Z
 Ribbon_Proofs
 Robbins-Conjecture
 Root_Balanced_Tree
 Routing
 Roy_Floyd_Warshall
 SATSolverVerification
 SDS_Impossibility
 SIFPL
 SIFUM_Type_Systems
 SPARCv8
 Safe_OCL
 Saturation_Framework
 Secondary_Sylow
 Security_Protocol_Refinement
 Selection_Heap_Sort
 SenSocialChoice
 Separata
 Separation_Algebra
 Separation_Logic_Imperative_HOL
 SequentInvertibility
 Shivers-CFA
 ShortestPath
 Show
 Sigma_Commit_Crypto
 Signature_Groebner
 Simpl
 Simple_Firewall
 Simplex
 Skew_Heap
 Skip_Lists
 Slicing
 Sliding_Window_Algorithm
 Smooth_Manifolds
 Sort_Encodings
 Source_Coding_Theorem
 Special_Function_Bounds
 Splay_Tree
 Sqrt_Babylonian
 Stable_Matching
 Statecharts
 Stellar_Quorums
 Stern_Brocot
 Stewart_Apollonius
 Stirling_Formula
 Stochastic_Matrices
 Stone_Algebras
 Stone_Kleene_Relation_Algebras
 Stone_Relation_Algebras
 Store_Buffer_Reduction
 Stream-Fusion
 Stream_Fusion_Code
 Strong_Security
 Sturm_Sequences
 Sturm_Tarski
 Stuttering_Equivalence
 Subresultants
 Subset_Boolean_Algebras
 SumSquares
 SuperCalc
 Surprise_Paradox
 Symmetric_Polynomials
 Szpilrajn
 TESL_Language
 TLA
 Tail_Recursive_Functions
 Tarskis_Geometry
 Taylor_Models
 Timed_Automata
 Topology
 TortoiseHare
 Transcendence_Series_Hancl_Rucki
 Transformer_Semantics
 Transition_Systems_and_Automata
 Transitive-Closure
 Transitive-Closure-II
 Treaps
 Tree-Automata
 Tree_Decomposition
 Triangle
 Trie
 Twelvefold_Way
 Tycon
 Types_Tableaus_and_Goedels_God
 Universal_Turing_Machine
 UPF
 UPF_Firewall
 UpDown_Scheme
 UTP
 Valuation
 VectorSpace
 VeriComp
 Verified-Prover
 VerifyThis2018
 VerifyThis2019
 Vickrey_Clarke_Groves
 VolpanoSmith
 WHATandWHERE_Security
 WebAssembly
 Weight_Balanced_Trees
 Well_Quasi_Orders
 Winding_Number_Eval
 WOOT_Strong_Eventual_Consistency
 Word_Lib
 WorkerWrapper
 XML
 Zeta_Function
 Zeta_3_Irrational
 ZFC_in_HOL
 pGCL