diff --git a/packages/preview/irif/0.0.1/LICENSE b/packages/preview/irif/0.0.1/LICENSE
new file mode 100644
index 0000000000..a0d8cb7181
--- /dev/null
+++ b/packages/preview/irif/0.0.1/LICENSE
@@ -0,0 +1,14 @@
+MIT No Attribution
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of
+this software and associated documentation files (the "Software"), to deal in
+the Software without restriction, including without limitation the rights to
+use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
+the Software, and to permit persons to whom the Software is furnished to do so.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
+FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
+COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
+IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
+CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
diff --git a/packages/preview/irif/0.0.1/README.md b/packages/preview/irif/0.0.1/README.md
new file mode 100644
index 0000000000..a62a211fb9
--- /dev/null
+++ b/packages/preview/irif/0.0.1/README.md
@@ -0,0 +1,98 @@
+# irif
+A [Typst](https://typst.app/) package for Numerical Methods to find roots, gradients and areas of functions.
+Graph plotting for integrals built on top of [CeTZ](https://github.com/johannes-wolf/cetz).
+
+Examples can be found [here as a .pdf](examples/examples.pdf) and [here as .typ](examples/examples.typ).
+
+
+```typ
+#plot-integral(f_x:x=>calc.pow(calc.e,x)+1/(x - 2),
+             x0:0.5,
+             x1:1.8, 
+             size-x:3,
+             size-y:2,
+             label-a:"0.5",
+             label-b:"1.8",
+             n-strips:8)
+
+#nm-table-integrate(f_x:x=>calc.pow(calc.e,x)+1/(x - 2),
+                    x0:0.5,
+                    x1:1.8,
+                    show-diffs:false,
+                    n-rows:10)
+
+```
+
+Functions names are as follows: `nm-[operation type]-[method]`
+Creates the following functions:
+```typ
+#nm-differentiate-forward()
+#nm-differentiate-central()
+
+#nm-integrate-midpoint()
+#nm-integrate-trapezium()
+#nm-integrate-simpsons()
+
+#nm-iterate-FPI()
+#nm-iterate-relaxed-FPI()
+#nm-iterate-newton-raphson()
+#nm-iterate-secant()
+#nm-iterate-false-position()
+#nm-iterate-bisection()
+
+#nm-table-integrate()
+#nm-table-differentiate()
+#nm-table-iterate()
+
+#plot-integral()
+```
+
+## nm-differentiation
+Differentiation tools take the following inputs/defaults:
+```typ
+f_x:x=>x*x, x0:1, h:1, accuracy:12
+```
+The function will approximate $f'(x_0)$ using either the Forward Difference or Central Difference method.
+
+## nm-integrate
+Integration tools take the following inputs/defaults:
+```typ
+f_x:x=>x*x, x0:0, x1:1, accuracy:12, n:1
+```
+These will approximate $\int_{x 0}^{x 1}(f_x(x))\text{d}x $
+
+Trapezium rule will return $T_n$, Midpoint rule will return $M_n$ and Simpsons rule will return $S_{2n}$
+
+## nm-iterate
+These functions provide the root-finding methods Bisection, False Position, Secant, Newton Raphson, Fixed Point Iteration (FPI) and Relaxed Fixed Point Iteration.
+
+All functions take the parameter `f_x`, except FPI and Relaxed FPI, which take the function `g_x`.
+
+Bisection and False Position will return the region in which the root lies, whereas the others will return their nearest approximation.
+
+The parameter `return-all` can be made true to see every step of the iteration, rather than just the last one requested.
+
+*NOTE*: The Newton-Raphson method utilises the Central Difference gradient with `h = 0.0000000001` calculated to 15 decimal places. Thus it is not a true Newton-Raphson approximation.
+
+## nm-table
+These tables are useful for seeing or demonstrating convergence. They have been designed with the A Level Further Maths OCR MEI B specification in mind.
+
+By default, they will include the changing variable ($n$ or $h$ typically), as well as the approximation they have reached.
+
+Optionally, the differences between the estimates can be shown.
+
+If the differences are being shown, then the ratios between the differences can also be shown, with a customizeable order for iterative methods.
+
+*NOTE*: For exam purposes, differences and ratios are caluclated from the *table values*, which are rounded, rather than the greater precision stored values.
+
+## plot-integral
+Useful for demonstrating different types of numerical integration methods. Uses CeTZ to plot.
+
+The number of strips is customizable, and the method should be chosen from:
+
+`integral,mid,left,right,trapezium`
+
+If `integral` is chosen, or any non-listed input is given for the `method` variable, the function defaults to simply highlighting the area to be found.
+
+*NOTE*: Labels aren't in mathematical format, simply New Computer Modern Math font.
+
diff --git a/packages/preview/irif/0.0.1/examples/examples.pdf b/packages/preview/irif/0.0.1/examples/examples.pdf
new file mode 100644
index 0000000000..22c4de658d
Binary files /dev/null and b/packages/preview/irif/0.0.1/examples/examples.pdf differ
diff --git a/packages/preview/irif/0.0.1/examples/examples.typ b/packages/preview/irif/0.0.1/examples/examples.typ
new file mode 100644
index 0000000000..47039dbbe2
--- /dev/null
+++ b/packages/preview/irif/0.0.1/examples/examples.typ
@@ -0,0 +1,167 @@
+#import "@preview/irif:0.0.1" : *
+#set page(margin: (top:2cm,x:1cm,bottom:1cm))
+== NM-Integrate
+#table(columns:5,inset:8pt,
+table.header([Target],$n$,[\#Midpoint],[\#Trapezium],[\#Simpsons]),
+table.cell(rowspan:4,[$pi approx $ #calc.round(calc.pi,digits:8)]),
+[1],[#nm-integrate-midpoint(f_x:x=>4 / (1 + x*x),x0:0,x1:1,n:1,accuracy:6)],
+  [#nm-integrate-trapezium(f_x:x=>4 / (1 + x*x),x0:0,x1:1,n:1,accuracy:6)],
+  [#nm-integrate-simpsons(f_x:x=>4 / (1 + x*x),x0:0,x1:1,n:1,accuracy:6)],
+[2],[#nm-integrate-midpoint(f_x:x=>4 / (1 + x*x),x0:0,x1:1,n:2,accuracy:6)],
+  [#nm-integrate-trapezium(f_x:x=>4 / (1 + x*x),x0:0,x1:1,n:2,accuracy:6)],
+  [#nm-integrate-simpsons(f_x:x=>4 / (1 + x*x),x0:0,x1:1,n:2,accuracy:6)],
+[4],[#nm-integrate-midpoint(f_x:x=>4 / (1 + x*x),x0:0,x1:1,n:4,accuracy:6)],
+  [#nm-integrate-trapezium(f_x:x=>4 / (1 + x*x),x0:0,x1:1,n:4,accuracy:6)],
+  [#nm-integrate-simpsons(f_x:x=>4 / (1 + x*x),x0:0,x1:1,n:4,accuracy:6)],
+[8],[#nm-integrate-midpoint(f_x:x=>4 / (1 + x*x),x0:0,x1:1,n:8,accuracy:6)],
+  [#nm-integrate-trapezium(f_x:x=>4 / (1 + x*x),x0:0,x1:1,n:8,accuracy:6)],
+  [#nm-integrate-simpsons(f_x:x=>4 / (1 + x*x),x0:0,x1:1,n:8,accuracy:6)],
+table.cell(rowspan:4,[$ln 2 approx $#calc.round(calc.ln(2),digits:8)]),
+[1],[#nm-integrate-midpoint(f_x:x=>1/x,x0:1,x1:2,n:1,accuracy:6)],
+  [#nm-integrate-trapezium(f_x:x=>1/x,x0:1,x1:2,n:1,accuracy:6)],
+  [#nm-integrate-simpsons(f_x:x=>1/x,x0:1,x1:2,n:1,accuracy:6)],
+[2],[#nm-integrate-midpoint(f_x:x=>1/x,x0:1,x1:2,n:2,accuracy:6)],
+  [#nm-integrate-trapezium(f_x:x=>1/x,x0:1,x1:2,n:2,accuracy:6)],
+  [#nm-integrate-simpsons(f_x:x=>1/x,x0:1,x1:2,n:2,accuracy:6)],
+[4],[#nm-integrate-midpoint(f_x:x=>1/x,x0:1,x1:2,n:4,accuracy:6)],
+  [#nm-integrate-trapezium(f_x:x=>1/x,x0:1,x1:2,n:4,accuracy:6)],
+  [#nm-integrate-simpsons(f_x:x=>1/x,x0:1,x1:2,n:4,accuracy:6)],
+[8],[#nm-integrate-midpoint(f_x:x=>1/x,x0:1,x1:2,n:8,accuracy:6)],
+  [#nm-integrate-trapezium(f_x:x=>1/x,x0:1,x1:2,n:8,accuracy:6)],
+  [#nm-integrate-simpsons(f_x:x=>1/x,x0:1,x1:2,n:8,accuracy:6)])
+
+== NM-Differentiate
+#table(columns:4,inset:8pt,
+table.header([Target],$h$,[\#Forward],[\#Central]),
+table.cell(rowspan:3,[$f(x)=cos x, quad f`(pi/2)=-1$]),
+$1$,[#nm-differentiate-forward(f_x:x=>calc.cos(x),x0:calc.pi/2,h:1,accuracy:6)],
+      [#nm-differentiate-central(f_x:x=>calc.cos(x),x0:calc.pi/2,h:1,accuracy:6)],
+$0.1$,[#nm-differentiate-forward(f_x:x=>calc.cos(x),x0:calc.pi/2,h:0.1,accuracy:6)],
+      [#nm-differentiate-central(f_x:x=>calc.cos(x),x0:calc.pi/2,h:0.1,accuracy:6)],
+$0.01$,[#nm-differentiate-forward(f_x:x=>calc.cos(x),x0:calc.pi/2,h:0.01,accuracy:6)],
+      [#nm-differentiate-central(f_x:x=>calc.cos(x),x0:calc.pi/2,h:0.01,accuracy:6)],
+table.cell(rowspan:3,[$f(x)=1/(1+ln x),quad f`(2) approx -0.174414$]),
+$1$,[#nm-differentiate-forward(f_x:x=>1 / (1 + calc.ln(x)),x0:2,h:1,accuracy:6)],
+  [#nm-differentiate-central(f_x:x=>1 / (1 + calc.ln(x)),x0:2,h:1,accuracy:6)],
+$0.1$,[#nm-differentiate-forward(f_x:x=>1 / (1 + calc.ln(x)),x0:2,h:0.1,accuracy:6)],
+  [#nm-differentiate-central(f_x:x=>1 / (1 + calc.ln(x)),x0:2,h:0.1,accuracy:6)],
+$0.01$,[#nm-differentiate-forward(f_x:x=>1 / (1 + calc.ln(x)),x0:2,h:0.01,accuracy:6)],
+  [#nm-differentiate-central(f_x:x=>1 / (1 + calc.ln(x)),x0:2,h:0.01,accuracy:6)])
+
+#pagebreak()
+== nm-iterate
+=== FPI
+#table(columns:4, inset:8pt,
+table.header([Target],[\#Fixed],[\#Relaxed $lambda=0.5$],[\#Relaxed $lambda=1.4$]),
+[$g(x)=2sin x + 2 cos x \ x approx -2.68075641016$],
+  [#nm-iterate-FPI(g_x:x=>2 * calc.sin(x) + 2 * calc.cos(x),accuracy:6,x0:-1.5,n-max:8)],
+  [#nm-iterate-relaxed-FPI(g_x:x=>2 * calc.sin(x) + 2 * calc.cos(x),accuracy:6,x0:-1.5,n-max:8,lambda:0.5)],
+  [#nm-iterate-relaxed-FPI(g_x:x=>2 * calc.sin(x) + 2 * calc.cos(x),accuracy:6,x0:-1.5,n-max:8,lambda:1.4)])
+All Values:
+#let FPI = nm-iterate-FPI(g_x:x=>2 * calc.sin(x) + 2 * calc.cos(x),accuracy:6,x0:-1.5,n-max:5,return-all:true)
+#let RFPI1 = nm-iterate-relaxed-FPI(g_x:x=>2 * calc.sin(x) + 2 * calc.cos(x),accuracy:6,x0:-1.5,n-max:5,lambda:0.5,return-all:true)
+#let RFPI2 = nm-iterate-relaxed-FPI(g_x:x=>2 * calc.sin(x) + 2 * calc.cos(x),accuracy:6,x0:-1.5,n-max:5,lambda:1.4,return-all:true)
+
+#let table-rows = ()
+#for i in range(FPI.len()) {
+  table-rows.push([#i])
+  table-rows.push([#FPI.at(i)])
+  table-rows.push([#RFPI1.at(i)])
+  table-rows.push([#RFPI2.at(i)])
+}
+
+  #table(columns:4,inset:5pt,
+  table.header($n$,[FPI $x_n$],[RFPI,$lambda=0.5$],[RFPI,$lambda=1.4$],
+  ..table-rows
+  )
+)
+=== Newton-Raphson
+#table(columns:2,inset:8pt,
+table.header([Target],[\#NRaphson]),
+[$f(x)=e^x-x-2 \ f(1.14619) approx 0 \ f(-1.84141) approx 9$],
+[#nm-iterate-newton-raphson(f_x:x=>calc.pow(calc.e,x) - x - 2,x0:1,accuracy:6) \ \
+  #nm-iterate-newton-raphson(f_x:x=>calc.pow(calc.e,x) - x - 2,x0:-1,accuracy:6)])
+All Values:
+#let NR1 = nm-iterate-newton-raphson(f_x:x=>calc.pow(calc.e,x) - x - 2,x0:1,accuracy:6,return-all:true)
+#let NR2 = nm-iterate-newton-raphson(f_x:x=>calc.pow(calc.e,x) - x - 2,x0:-1,accuracy:6,return-all:true)
+#let table-rows = ()
+#for i in range(NR1.len()) {
+  table-rows.push([#i])
+  table-rows.push([#NR1.at(i)])
+  table-rows.push([#NR2.at(i)])
+}
+#table(columns:3,inset:5pt,
+table.header([$n$],[NR:$x_0 = 1$],[NR:$x_0 = -1$]),
+..table-rows
+)
+
+=== Secant, False Position, Bisection for $sin(x)=0$
+#table(columns:(0.2fr,1fr,1fr,1fr),inset:5pt,
+table.header($x$,[Secant],[False Position],[Bisection]),
+text(size:8pt)[$[-1,1]$],
+  [#nm-iterate-secant(x0:-1,x1:1,return-all:true,accuracy:8)],
+  [#nm-iterate-false-position(x0:-1,x1:1,return-all:true,accuracy:8)],
+  [#nm-iterate-bisection(x0:-1,x1:1,return-all:true,accuracy:8)],
+text(size:8pt)[$[0.5,1]$],
+[#nm-iterate-secant(x0:0.5,x1:1,return-all:true,accuracy:8)],
+  [ False Position should panic.
+  //#nm-iterate-false-position(x0:0.5,x1:1,return-all:true,accuracy:8)
+],
+  [ Bisection should panic.
+  //#nm-iterate-bisection(x0:0.5,x1:1,return-all:true,accuracy:8)
+],
+text(size:8pt)[$[-1.5,1]$],
+  [#nm-iterate-secant(x0:-1.5,x1:1,return-all:true,accuracy:8)],
+  [#nm-iterate-false-position(x0:-1.5,x1:1,return-all:true,accuracy:8)],
+  [#nm-iterate-bisection(x0:-1.5,x1:1,return-all:true,accuracy:8)],
+)
+
+#pagebreak()
+== nm-plot-integral
+=== Types: Left, Right, Midpoint, Trapezium, Integral, 'asfsf'
+Alternatively including 'points'
+#grid(columns:3, inset:10pt,
+plot-integral(size-x:3,size-y:2,method:"left",show-points:true), 
+plot-integral(size-x:3,size-y:2,method:"right"),
+plot-integral(size-x:3,size-y:2,method:"mid",show-points:true), 
+plot-integral(size-x:3,size-y:2,method:"trapezium"),
+plot-integral(size-x:3,size-y:2,method:"integral",show-points:true), 
+plot-integral(size-x:3,size-y:2,method:"asfsf"), 
+)
+=== Limits: inc_x0, no x0
+#grid(columns:3,inset:10pt,
+plot-integral(size-x:3,size-y:2,method:"integral",x0:2,x1:3), 
+plot-integral(size-x:3,size-y:2,method:"integral",x0:2,x1:3,inc-0:false),
+)
+=== Strips: 1, 2, 4, 8, 16, 32
+#grid(columns:3, inset:10pt,
+plot-integral(size-x:3,size-y:2,n-strips:1), 
+plot-integral(size-x:3,size-y:2,n-strips:2), 
+plot-integral(size-x:3,size-y:2,n-strips:4), 
+plot-integral(size-x:3,size-y:2,n-strips:8), 
+plot-integral(size-x:3,size-y:2,n-strips:16), 
+plot-integral(size-x:3,size-y:2,n-strips:32)
+)
+
+#pagebreak()
+== NM-Table
+=== Integral Table
+#grid(columns:3,inset: 5pt,
+nm-table-integrate(),
+nm-table-integrate(method:"Trapezium"),
+nm-table-integrate(method:"Simpsons",n-rows:4)
+)
+
+=== Differential Table
+#grid(columns:2,inset: 5pt,
+nm-table-differentiate(),
+nm-table-differentiate(method:"FD"),
+)
+
+=== Iteration Table
+#grid(columns:2,inset: 5pt,
+[FPI: #nm-table-iterate(method:"FPI",ratio-order:1)],
+[RFPI: #nm-table-iterate(method:"RFPI",ratio-order:1)],
+[NRaphson: #nm-table-iterate(method:"Newton-Raphson",ratio-order:2)],
+[Secant: #nm-table-iterate(method:"Secant",ratio-order:1.6)]
+)
\ No newline at end of file
diff --git a/packages/preview/irif/0.0.1/irif.typ b/packages/preview/irif/0.0.1/irif.typ
new file mode 100644
index 0000000000..ae91f770dc
--- /dev/null
+++ b/packages/preview/irif/0.0.1/irif.typ
@@ -0,0 +1,468 @@
+#import "@preview/cetz:0.4.2"
+#import "@preview/cetz-plot:0.1.3"
+
+#let plot-integral(f_x: x => x * x,
+                  x0: 0.2,
+                  x1: 6,
+                  inc-0: true,
+                  shade-color: rgb(100,100,255,50),
+                  method:"mid",
+                  n-strips:5,
+                  show-points:false,
+                  show-labels:true,
+                  label-a: "a",
+                  label-b:"b",
+                  size-x:6,size-y:6
+) = {
+  // Mistaken plot type:
+  if (not ("integral","mid","left","right","trapezium").contains(method)) {method = "integral"}
+  // Fix incorrect bounds:
+  if (x0 > x1) {(x0,x1) = (x1,x0)}
+  else if (x0==x1){x1 = x0 + 1}
+  //Create variables
+  let min-x = 0
+  let max-x = 0
+  let min-y = 0
+  let max-y = 0
+  let h = (x1 - x0) / n-strips // Strip width, if used.
+  if (inc-0) { //Actually calculate X-min/max
+    let width = x1
+    max-x = width * 1.1
+    min-x = - width * 0.05 
+  } else {
+    let width = x1 - x0
+    max-x = x1 + 0.05 * width
+    min-x = x0 - 0.05 * width
+  }
+  //Largest found y-value
+  let f-max = calc.max(..range(0,1001).map(x=>f_x(x0 + x*(x1 - x0)/1000)))
+  //values for plot
+  let min-y = -0.05 * f-max
+  let max-y = 1.1 * f-max
+  
+  
+  cetz.canvas({
+    import cetz.draw: *
+    import cetz-plot: *
+    set-style(stroke:2pt)
+    plot.plot(size:(size-x,size-y),x-tick-step: none,y-tick-step: none, axis-style:"school-book",
+    x-min:min-x,x-max:max-x,y-min:min-y,y-max:max-y,
+    {
+
+      plot.add(domain:(min-x,max-x), f_x,style:(stroke:black))
+      
+      if (method=="integral") {
+        plot.add(domain:(x0,x1),f_x,fill-type:"axis",fill:true,style:(fill:shade-color,stroke:3pt))
+        plot.add-vline(x0,min:0,max:f_x(x0),style:(stroke:1pt+black))
+        plot.add-vline(x1,min:0,max:f_x(x1),style:(stroke:1pt+black))
+        if (show-points) {
+            plot.add(((x0,f_x(x0)),(x0,f_x(x0))),mark:"*",mark-style:(stroke:black))
+            plot.add(((x1,f_x(x1)),(x1,f_x(x1))),mark:"*",mark-style:(stroke:black))
+          }
+      } else if (method=="left"){
+        let n = 0
+        let x = x0
+        while n < n-strips {
+          plot.add(domain:(x,x+h),p=>f_x(x),style:(fill:shade-color,stroke:1pt+black),fill:true)
+          plot.add-vline(x,min:0,max:f_x(x),style:(stroke:1pt+black))
+          plot.add-vline(x+h,min:0,max:f_x(x),style:(stroke:1pt+black))
+          if (show-points) {
+            plot.add(((x,f_x(x)),(x,f_x(x))),mark:"*",mark-style:(stroke:black))
+          }
+          x += h
+          n += 1
+        }
+      } else if (method=="right"){
+       let n = 0
+        let x = x0
+        while n < n-strips {
+          plot.add(domain:(x,x+h),p=>f_x(x+h),style:(fill:shade-color,stroke:1pt+black),fill:true)
+          plot.add-vline(x,min:0,max:f_x(x+h),style:(stroke:1pt+black))
+          plot.add-vline(x+h,min:0,max:f_x(x+h),style:(stroke:1pt+black))
+          if (show-points) {
+            plot.add(((x+h,f_x(x+h)),(x+h,f_x(x+h))),mark:"*",mark-style:(stroke:black))
+          }
+          x += h
+          n += 1
+        }
+      } else if (method=="mid") {
+        let n = 0
+        let x = x0
+        while n < n-strips {
+          plot.add(domain:(x,x+h),p=>f_x(x+h/2),style:(fill:shade-color,stroke:1pt+black),fill:true)
+          plot.add-vline(x,min:0,max:f_x(x+h/2),style:(stroke:1pt+black))
+          plot.add-vline(x+h,min:0,max:f_x(x+h/2),style:(stroke:1pt+black))
+          if (show-points) {
+            plot.add(((x+h/2,f_x(x+h/2)),(x+h/2,f_x(x+h/2))),mark:"*",mark-style:(stroke:black))
+          }
+          x += h
+          n += 1
+        }
+      } else if (method=="trapezium"){
+        let n = 0
+        let x = x0
+        while n < n-strips {
+          plot.add(domain:(x,x+h),n=>(n - x)*(f_x(x+h)-f_x(x))/(h)+f_x(x),style:(fill:shade-color,stroke:1pt+black),fill:true)
+          plot.add-vline(x,min:0,max:f_x(x),style:(stroke:1pt+black))
+          plot.add-vline(x+h,min:0,max:f_x(x+h),style:(stroke:1pt+black))
+          if (show-points) {
+            plot.add(((x,f_x(x)),(x,f_x(x))),mark:"*",mark-style:(stroke:black))
+          }
+          x += h
+          n += 1
+        }
+        if (show-points) {plot.add(((x,f_x(x)),(x,f_x(x))),mark:"*",mark-style:(stroke:black))}
+      }
+      if (show-labels){
+        plot.annotate({
+          content((x0,2*min-y),text(font:"New Computer Modern Math",label-a))
+          content((x1,2*min-y),text(font:"New Computer Modern Math",label-b))
+        })
+      }
+    })
+  })  
+}
+
+#let nm-integrate-midpoint(f_x:x=>x*x,
+x0:0,x1:1,accuracy:12,n:1) = {
+  let h = (x1 - x0) / n
+  let total = 0
+  for i in range(n) {
+    total += h * f_x(h/2 + x0 + h * i)
+  }
+  return calc.round(total,digits:accuracy)
+}
+#let nm-integrate-trapezium(f_x:x=>x*x,
+x0:0,x1:1,accuracy:12,n:1) = {
+  let h = (x1 - x0) / n
+  let total = 0.5 * (f_x(x0)+f_x(x1))
+  let strip = 1
+  while strip < n {
+    total += f_x(x0 + strip * h)
+    strip += 1
+  }
+  return calc.round(h*total,digits:accuracy)
+}
+#let nm-integrate-simpsons(f_x:x=>x*x,
+x0:0,x1:1,accuracy:12,n:1) = {
+  let T_n = nm-integrate-trapezium(f_x:f_x,x0:x0,x1:x1,accuracy:accuracy+5,n:n)
+  let T_2n = nm-integrate-trapezium(f_x:f_x,x0:x0,x1:x1,accuracy:accuracy+5,n:2*n)
+
+  let total = (4 * T_2n - T_n) / 3
+  
+  return calc.round(total,digits:accuracy)
+}
+
+#let nm-differentiate-forward(f_x:x=>x*x,x0:1,h:1,accuracy:12) = {
+  let d = (f_x(x0+h) - f_x(x0)) / h
+  return calc.round(d,digits:accuracy)
+}
+#let nm-differentiate-central(f_x:x=>x*x,x0:1,h:1,accuracy:12) = {
+  let d = (f_x(x0 + h/2) - f_x(x0 - h/2)) / h
+  return calc.round(d,digits:accuracy)
+}
+
+
+
+#let nm-iterate-relaxed-FPI(g_x:calc.cos,x0:1,lambda:0.5,n-max:5,accuracy:12,return-all:false) = {
+  let n = 0
+
+  if not return-all {
+    let x = x0
+    while n < n-max {
+      x = lambda * g_x(x) + (1 - lambda) * x
+      n += 1
+    }
+    return calc.round(x,digits:accuracy)
+  }
+  
+  let x = ()
+  x.push(x0)
+  while n < n-max {
+    x.push(lambda * g_x(x.last()) + (1 - lambda)*x.last())
+    n += 1
+  }
+  return x.map(n=>calc.round(n,digits:accuracy))
+  
+}
+#let nm-iterate-FPI(g_x:calc.cos,x0:1,n-max:5,accuracy:12,return-all:false) = {
+  let result = nm-iterate-relaxed-FPI(g_x:g_x,x0:x0,n-max:n-max,lambda:1,accuracy:accuracy,return-all:return-all)
+  return result
+}
+#let nm-iterate-newton-raphson(f_x:calc.cos,x0:1,n-max:5,accuracy:12,return-all:false) = {
+  let n = 0
+
+  if not return-all {
+    let x = x0
+    while n < n-max {
+      x = x - (f_x(x)) / (nm-differentiate-central(f_x:f_x,
+                                            x0:x,
+                                            h:0.0000001,
+                                            accuracy:15))
+      n += 1
+    }
+    return calc.round(x,digits:accuracy)
+  }
+  let x = ()
+  x.push(x0)
+  while n < n-max {
+    let xn = x.last()
+    x.push(xn - (f_x(xn)) / (nm-differentiate-central(f_x:f_x,
+                                            x0:xn,
+                                            h:0.0000000001,
+                                            accuracy:15)))
+    n += 1
+  }
+  return x.map(n=>calc.round(n,digits:accuracy))
+}
+
+#let nm-iterate-secant(f_x:calc.sin,x0:1,x1:0.5,n-max: 5,accuracy:5,return-all:false) = {
+  let n = 0
+  if not return-all {
+    let (a,b) = (x0,x1)
+    while n < n-max {
+      let next = (a * f_x(b) - b * f_x(a))/(f_x(b) - f_x(a))
+      a = b
+      b = next
+      n += 1
+    }
+    return calc.round(b,digits:accuracy)
+  }
+  let xn = ()
+  xn.push(x0)
+  xn.push(x1)
+  while xn.len() < n-max {
+    let (a,b) = (xn.at(xn.len()-2),xn.last())
+    if (f_x(b) - f_x(a))==0 {
+      xn.push(0)
+    } else {
+    let next = (a * f_x(b) - b * f_x(a))/(f_x(b) - f_x(a))
+    xn.push(next)
+    }
+  }
+  return xn.map(n=>calc.round(n,digits:accuracy))
+}
+
+#let nm-iterate-false-position(f_x:calc.sin,x0:1,x1:0.5,n-max: 5,accuracy:5,return-all:false) = {
+  if f_x(x0) * f_x(x1) > 0 {panic("Error, no change of sign.")}
+  let n = 0
+  if not return-all {
+    let (a,b) = (x0,x1)
+    while n < n-max {
+      if (f_x(a) * f_x(b) >=0) {panic("Error, no change of sign.")}
+      let next = (a * f_x(b) - b * f_x(a))/(f_x(b) - f_x(a))
+      if f_x(a)*f_x(next) > 0 { //Same sign, so replace b
+        b = next
+      } else {a = next}
+      n += 1
+    }
+    return (calc.round(a, digits: accuracy), 
+            calc.round(b, digits: accuracy))
+  }
+
+  let xn = ()
+  xn.push((x0,x1))
+  while xn.len() < n-max {
+    let (a,b) = xn.last()
+    if (f_x(a) * f_x(b) >=0) {panic("Error, no change of sign.")}
+    let next = (a * f_x(b) - b * f_x(a))/(f_x(b) - f_x(a))
+    if f_x(a)*f_x(next) > 0 { //Same sign, so replace b
+        xn.push((a,next))
+      } else {xn.push((next,b))}
+    n += 1
+  }
+  return xn.map(n=>(calc.round(n.at(0), digits: accuracy),
+                        calc.round(n.at(1), digits: accuracy)))
+}
+#let nm-iterate-bisection(f_x:calc.sin,x0:1,x1:0.5,n-max: 5,accuracy:5,return-all:false) = {
+  if f_x(x0) * f_x(x1) > 0 {panic("Error, no change of sign.")}
+  let n = 0
+  if not return-all {
+    let (a,b) = (x0,x1)
+    while n < n-max {
+      if (f_x(a) * f_x(b) >=0) {panic("Error, no change of sign.")}
+      let next = (a+b)/2
+      if (f_x(next)==0) {
+      xn.push((a,next))
+      break
+    }
+      if f_x(a)*f_x(next) > 0 { //Same sign, so replace b
+        b = next
+      } else {a = next}
+      n += 1
+    }
+    return (calc.round(a, digits: accuracy), 
+            calc.round(b, digits: accuracy))
+  }
+
+  let xn = ()
+  xn.push((x0,x1))
+  while xn.len() < n-max {
+    let (a,b) = xn.last()
+    if (f_x(a) * f_x(b) >=0) {panic("Error, no change of sign.")}
+    let next = (a+b)/2
+    if (f_x(next)==0) {
+      xn.push((a,next))
+      break
+    }
+    if f_x(a)*f_x(next) > 0 { //Same sign, so replace b
+        xn.push((a,next))
+      } else {xn.push((next,b))}
+    n += 1
+  }
+  return xn.map(n=>(calc.round(n.at(0), digits: accuracy),
+                        calc.round(n.at(1), digits: accuracy)))
+}
+
+#let nm-table-integrate(f_x:calc.sin,x0:1,x1:2,start-strips:1,strip-ratio:2,method:"Midpoint", n-rows:5,show-diffs:true,show-ratio:true,accuracy:5) = {
+  //method must be from the array. If not, do mid.
+  if not ("Midpoint","Trapezium","Simpsons").contains(method) {
+    method = "Midpoint"
+  }
+  let strip-array = ()
+  strip-array.push(start-strips)
+  while strip-array.len() < n-rows {
+    strip-array.push(strip-array.last() * strip-ratio)
+  }
+
+  let integrals = ()
+  let header = $M_n$
+  if method == "Midpoint" {
+    integrals = strip-array.map(n=>nm-integrate-midpoint(f_x:f_x,x0:x0,x1:x1,n:n,accuracy:accuracy))
+  } else if method=="Trapezium" {
+    integrals = strip-array.map(n=>nm-integrate-trapezium(f_x:f_x,x0:x0,x1:x1,n:n,accuracy:accuracy))
+    header = $T_n$
+  } else {
+    integrals = strip-array.map(n=>nm-integrate-simpsons(f_x:f_x,x0:x0,x1:x1,n:n,accuracy:accuracy))
+    header = $S_(2n)$
+  }
+  let diffs = ()
+  diffs.push([])
+  for i in range(integrals.len()-1) {
+    diffs.push(calc.round(integrals.at(i+1) - integrals.at(i),digits:accuracy))
+  }
+  let ratios = ([],[])
+  for i in range(diffs.len()-2) {
+    ratios.push(calc.round(diffs.at(i+2) / diffs.at(i+1),digits:accuracy))
+  }
+
+  let table-rows = () 
+  for i in range(strip-array.len()) {
+    table-rows.push([#strip-array.at(i)])
+    table-rows.push([#integrals.at(i)])
+    if show-diffs {
+      table-rows.push([#diffs.at(i)])
+      if show-ratio {table-rows.push([#ratios.at(i)])}
+    }
+  }
+  let headers = ($n$,[#header])
+  if show-diffs {
+    headers.push([Difference])
+    if show-ratio {headers.push([Ratio])}
+  }
+  return table(columns:headers.len(),
+              table.header(..headers),
+              ..table-rows)
+}
+
+#let nm-table-differentiate(f_x:x=>calc.ln(x),x0:2,start-h:1,h-ratio:2,n-rows:5,accuracy:5,method:"CD",show-diffs:true,show-ratio:true) = {
+  //Default to Central Difference if bad name
+  if not ("CD","FD").contains(method) {
+    method = "CD"
+  }
+
+  let differentials = ()
+
+  let h-array = ()
+  h-array.push(start-h)
+  while h-array.len() < n-rows {
+    h-array.push(calc.round(h-array.last() / h-ratio,digits:8)) 
+    //Rounds due to floating point division
+  }
+  
+  let differentials = h-array.map(h=>nm-differentiate-central(f_x:f_x,x0:x0,h:h,accuracy:accuracy))
+  if method == "FD" {
+    differentials = h-array.map(h=>nm-differentiate-forward(f_x:f_x,x0:x0,h:h,accuracy:accuracy))
+  }
+
+  let diffs = ()
+  diffs.push([])
+  for i in range(differentials.len()-1) {
+    diffs.push(calc.round(differentials.at(i+1) - differentials.at(i),digits:accuracy))
+  }
+  let ratios = ([],[])
+  for i in range(diffs.len()-2) {
+    ratios.push(calc.round(diffs.at(i+2) / diffs.at(i+1),digits:accuracy))
+  }
+
+  let table-rows = () 
+  for i in range(h-array.len()) {
+    table-rows.push([#h-array.at(i)])
+    table-rows.push([#differentials.at(i)])
+    if show-diffs {
+      table-rows.push([#diffs.at(i)])
+      if show-ratio {table-rows.push([#ratios.at(i)])}
+    }
+  }
+  let headers = ($h$,$f'(x)$)
+  if show-diffs {
+    headers.push([Difference])
+    if show-ratio {headers.push([Ratio])}
+  }
+  return table(columns:headers.len(),
+              table.header(..headers),
+              ..table-rows)
+  
+}
+
+#let nm-table-iterate(f_x:calc.sin,x0:1,x1:0.5,lambda:0.5,n-rows: 5,accuracy:5,ratio-order:2,method:"FPI",show-diffs:true,show-ratio:true) = {
+  //Default to NRaphson
+  if not ("FPI","RFPI","Secant","Newton-Raphson").contains(method) {
+    method = "Newton-Raphson"
+  }
+
+  let n-array = range(n-rows).map(n=>n+1) //Start at 1.
+
+  let headers = ($n$,$x_n$)
+  if show-diffs {
+    headers.push([Difference])
+    if show-ratio {
+      headers.push([Ratio])
+    }
+  }
+  let x-array = ()
+  if method=="FPI" {
+    x-array = nm-iterate-FPI(g_x:f_x,x0:x0,n-max:n-rows,return-all:true,accuracy:accuracy)
+  } else if method=="RFPI" {
+    x-array = nm-iterate-relaxed-FPI(g_x:f_x,x0:x0,n-max:n-rows,return-all:true,accuracy:accuracy,lambda:lambda)
+  }else if method == "Secant" {
+    x-array = nm-iterate-secant(f_x:f_x,x0:x0,x1:x1,n-max:n-rows,return-all:true,accuracy:accuracy)
+  } else {
+    x-array = nm-iterate-newton-raphson(f_x:f_x,x0:x0,n-max:n-rows,return-all:true,accuracy:accuracy)
+  }
+  let diffs = ()
+  diffs.push([])
+  for i in range(x-array.len()-1) {
+    diffs.push(calc.round(x-array.at(i+1) - x-array.at(i),digits:accuracy))
+  }
+  let ratios = ([],[])
+  for i in range(diffs.len()-2) {
+    ratios.push(calc.round(calc.abs(diffs.at(i+2)) / calc.pow(calc.abs(diffs.at(i+1)),ratio-order),digits:accuracy))
+  }
+  
+  let table-rows = ()
+  for i in range(n-array.len()) {
+    table-rows.push([#n-array.at(i)])
+    table-rows.push([#x-array.at(i)])
+    if show-diffs {
+      table-rows.push([#diffs.at(i)])
+      if show-ratio {table-rows.push([#ratios.at(i)])}
+    }
+  }
+
+  return table(columns:headers.len(),
+              table.header(..headers),
+              ..table-rows)
+
+  
+}
\ No newline at end of file
diff --git a/packages/preview/irif/0.0.1/typst.toml b/packages/preview/irif/0.0.1/typst.toml
new file mode 100644
index 0000000000..2eeea98d23
--- /dev/null
+++ b/packages/preview/irif/0.0.1/typst.toml
@@ -0,0 +1,12 @@
+[package]
+name        = "irif"
+version     = "0.0.1"
+description = "Numerical methods for integration, differentiation and root finding"
+license     = "MIT-0"
+authors     = ["TGGrace"]
+entrypoint  = "irif.typ"
+exclude     = [
+    ".git*",
+    "examples"
+]
+